Version:	4.21.0
Date:	2024-11-10
Title:	Structural Equation Modeling and Twin Modeling in R
Maintainer:	Timothy C. Bates <timothy.c.bates@gmail.com>
License:	GPL-3
Language:	en-US
Encoding:	UTF-8
URL:	https://github.com/tbates/umx#readme
Description:	Quickly create, run, and report structural equation models, and twin models. See '?umx' for help, and umx_open_CRAN_page("umx") for NEWS. Timothy C. Bates, Michael C. Neale, Hermine H. Maes, (2019). umx: A library for Structural Equation and Twin Modelling in R. Twin Research and Human Genetics, 22, 27-41. <doi:10.1017/thg.2019.2>.
Depends:	R (≥ 4.1.0), OpenMx (≥ 2.19.0),
Imports:	cowplot, DiagrammeR, ggplot2, kableExtra, knitr, lavaan, MASS, Matrix, methods, MuMIn, mvtnorm, nlme, polycor, R2HTML, RCurl, scales, utils, xtable, gert, hrbrthemes, openxlsx, paran, psych, foreign, psychTools, pwr, rmarkdown
Suggests:	cocor, devtools, rhub, spelling, testthat, parallel, GPArotation
Enhances:	DiagrammeRsvg, rsvg
BugReports:	https://github.com/tbates/umx/issues
LazyData:	true
RoxygenNote:	7.3.2
NeedsCompilation:	no
Packaged:	2024-11-15 01:16:35 UTC; timothybates
Author:	Timothy C. Bates [aut, cre], Nathan Gillespie [wit], Hermine Maes [ctb], Michael C. Neale [ctb], Joshua N. Pritikin [ctb], Luis De Araujo [ctb], Brenton Wiernik [ctb], Michael Zakharin [wit]
Repository:	CRAN
Date/Publication:	2024-11-15 03:30:02 UTC

Functions for Structural Equation Modeling in OpenMx

Description

umx allows you to more easily build, run, modify, and report structural models, building on the OpenMx package. All core functions are organized into families, so they are easier to find (so if you know a function similar to what you are looking for, look at other members of its "family" at the bottom of its help file.

Please cite as: Bates, T. C., Neale, M. C., & Maes, H. H. (2019). umx: A library for Structural Equation and Twin Modelling in R. Twin Research and Human Genetics, 22, 27-41. doi:10.1017/thg.2019.2.

All the functions have full-featured and well commented examples, some even have figures, so use the help, even if you think it won't help :-) Have a look, for example at umxRAM()

Check out NEWS about new features at news(package = "umx")

Details

Introductory working examples are below. You can run all demos with demo(umx) When I have a vignette, it will be: vignette("umx", package = "umx")

There is a helpful blog at https://tbates.github.io

(Only) if you want the bleeding-edge version:

devtools::install_github("tbates/umx")

Author(s)

Maintainer: Timothy C. Bates timothy.c.bates@gmail.com (ORCID)

Other contributors:

• Nathan Gillespie nathan.gillespie@vcuhealth.org [witness]

• Hermine Maes hmeaes@vcu.edu [contributor]

• Michael C. Neale neale@vcu.edu [contributor]

• Joshua N. Pritikin jpritikin@pobox.com [contributor]

• Luis De Araujo ldearaujo@unimelb.edu.au [contributor]

• Brenton Wiernik wiernik@umn.edu [contributor]

• Michael Zakharin s1775682@sms.ed.ac.uk [witness]

References

• Bates, T. C., Neale, M. C., & Maes, H. H. (2019). umx: A library for Structural Equation and Twin Modelling in R. Twin Research and Human Genetics, 22, 27-41. doi:10.1017/thg.2019.2, https://github.com/tbates/umx, tutorial: https://tbates.github.io

See Also

Useful links:

• https://github.com/tbates/umx#readme

• Report bugs at https://github.com/tbates/umx/issues

Other Core Model Building Functions: umxMatrix(), umxModify(), umxPath(), umxRAM(), umxSuperModel()

Other Model Summary and Comparison: umxCompare(), umxEquate(), umxMI(), umxReduce(), umxSetParameters(), umxSummary()

Other Reporting Functions: umxAPA(), umxFactorScores(), umxGetLatents(), umxGetManifests(), umxGetModel(), umxGetParameters(), umxParameters(), umx_aggregate(), umx_time()

Other Super-easy helpers: umxEFA(), umxTwoStage()

Other Twin Modeling Functions: power.ACE.test(), umxACE(), umxACEcov(), umxACEv(), umxCP(), umxDiffMZ(), umxDiscTwin(), umxDoC(), umxDoCp(), umxGxE(), umxGxE_window(), umxGxEbiv(), umxIP(), umxMRDoC(), umxReduce(), umxReduceACE(), umxReduceGxE(), umxRotate.MxModelCP(), umxSexLim(), umxSimplex(), umxSummarizeTwinData(), umxSummaryACE(), umxSummaryACEv(), umxSummaryDoC(), umxSummaryGxEbiv(), umxSummarySexLim(), umxSummarySimplex(), umxTwinMaker()

Other Twin Data functions: umx_long2wide(), umx_make_TwinData(), umx_make_twin_data_nice(), umx_residualize(), umx_scale_wide_twin_data(), umx_wide2long()

Other Miscellaneous Stats Functions: FishersMethod(), SE_from_p(), geometric_mean(), harmonic_mean(), oddsratio(), reliability(), umxCov2cor(), umxHetCor(), umxParan(), umxWeightedAIC(), umx_apply(), umx_cor(), umx_means(), umx_r_test(), umx_round(), umx_scale(), umx_var()

Other Teaching and testing Functions: tmx_genotypic_effect(), tmx_is.identified()

Other Get and set: umx_get_checkpoint(), umx_get_options(), umx_set_auto_plot(), umx_set_auto_run(), umx_set_checkpoint(), umx_set_condensed_slots(), umx_set_cores(), umx_set_data_variance_check(), umx_set_dollar_symbol(), umx_set_optimization_options(), umx_set_optimizer(), umx_set_plot_file_suffix(), umx_set_plot_format(), umx_set_plot_use_hrbrthemes(), umx_set_separator(), umx_set_silent(), umx_set_table_format()

Other Check or test: umx_check_names(), umx_is_class(), umx_is_endogenous(), umx_is_exogenous(), umx_is_numeric(), umx_is_ordered()

Other Plotting functions: ggAddR(), plot.MxLISRELModel(), plot.MxModel(), plot.MxModelTwinMaker(), umxPlot(), umxPlotACE(), umxPlotACEcov(), umxPlotACEv(), umxPlotCP(), umxPlotDoC(), umxPlotFun(), umxPlotGxE(), umxPlotGxEbiv(), umxPlotIP(), umxPlotSexLim(), umxPlotSimplex()

Other Data Functions: noNAs(), prolific_anonymize(), prolific_check_ID(), prolific_read_demog(), umxFactor(), umxHetCor(), umx_as_numeric(), umx_cont_2_quantiles(), umx_lower2full(), umx_make_MR_data(), umx_make_TwinData(), umx_make_fake_data(), umx_make_raw_from_cov(), umx_merge_randomized_columns(), umx_polychoric(), umx_polypairwise(), umx_polytriowise(), umx_read_lower(), umx_rename(), umx_reorder(), umx_score_scale(), umx_select_valid(), umx_stack(), umx_strings2numeric()

Other File Functions: dl_from_dropbox(), umx_file_load_pseudo(), umx_make_sql_from_excel(), umx_move_file(), umx_open(), umx_rename_file(), umx_write_to_clipboard()

Other String Functions: umx_explode(), umx_explode_twin_names(), umx_grep(), umx_names(), umx_paste_names(), umx_rot(), umx_str_chars(), umx_str_from_object(), umx_trim()

Other Miscellaneous Utility Functions: install.OpenMx(), libs(), qm(), umxLav2RAM(), umxModelNames(), umxRAM2Lav(), umxVersion(), umx_array_shift(), umx_find_object(), umx_lower.tri(), umx_msg(), umx_open_CRAN_page(), umx_pad(), umx_print()

Other datasets: Fischbein_wt, GFF, docData, iqdat, us_skinfold_data

Other Advanced Model Building Functions: umxAlgebra(), umxFixAll(), umxJiggle(), umxRun(), umxThresholdMatrix(), umxUnexplainedCausalNexus(), xmuLabel(), xmuValues()

Examples

## Not run: 
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 added informative labels, created starting values, 
# Ran your model (if autoRun is on), and displayed a brief summary
# including a comparison if you modified a model...!

# umxSummary generates journal-ready fit information.
# We can choose std=T for standardized parameters and can also
# filter out some types of parameter (e.g. means or residuals)

umxSummary(m1, std = TRUE, residuals=FALSE)

# parameters() flexibly retrieves model coefficients. 
# For example just G-loadings greater than |.3| and rounded to 2-digits.
parameters(m1, thresh="above", b=.3, pattern = "G_to.*", digits = 2)

# (The built-in coef works as for lm etc.)
coef(m1)

# ==================
# = Model updating =
# ==================
# umxModify modifies, renames, re-runs, and compares a model
# Can we set the loading of x1 on G to zero? (nope...)
m2 = umxModify(m1, "G_to_x1", name = "no_effect_of_g_on_X1", comparison = TRUE)

# note1: umxSetParameters can do this with some additional flexibility
# note2 "comparison = TRUE" above is the same as calling 
# umxCompare, like this
umxCompare(m1, m2)


# ========================
# = Confidence intervals =
# ========================

# umxSummary() will show these, but you can also use the confint() function
confint(m1) # OpenMx's SE-based confidence intervals


# umxConfint formats everything you need nicely, and allows adding CIs (with parm=)
umxConfint(m1, parm = 'all', run = TRUE) # likelihood-based CIs

# And make a Figure and open in browser
plot(m1, std = TRUE)

# If you just want the .dot code returned set file = NA
plot(m1, std = TRUE, file = NA)

## End(Not run)

Weight data across time.

Description

A dataframe containing correlations of weight for 66 females measured 6 times at 6-month intervals.

Usage

data(Fischbein_wt)

Format

A 6*6 correlation matrix based on n = 66 female subjects.

Details

• Weight1: Weight at time 1 (t0)

• Weight2: Weight at time 2 (t0 + 6 months)

• Weight3: Weight at time 3 (t0 + 12 months)

• Weight4: Weight at time 4 (t0 + 18 months)

• Weight5: Weight at time 5 (t0 + 24 months)

• Weight6: Weight at time 6 (t0 + 32 months)

Created as follows:

Fischbein_wt = umx_read_lower(file = "", diag = TRUE, names = paste0("Weight", 1:6), ensurePD= TRUE)
1.000
0.985	1.000
0.968	0.981	1.000
0.957	0.970	0.985	1.000
0.932	0.940	0.964	0.975	1.000
0.890	0.897	0.927	0.949	0.973	1.000

References

Fischbein, S. (1977). Intra-pair similarity in physical growth of monozygotic and of dizygotic twins during puberty. Annals of Human Biology, 4. 417-430. doi:10.1080/03014467700002401

See Also

Other datasets: GFF, docData, iqdat, umx, us_skinfold_data

Examples

## Not run: 
data(Fischbein_wt) # load the data
str(Fischbein_wt) # data.frame
as.matrix(Fischbein_wt) # convert to matrix

## End(Not run)

Fishers Method of combining p-values.

Description

FishersMethod implements R.A. Fisher's (1925) method for creating a meta-analytic p-value by combining a set of p-values from tests of the same hypothesis in independent samples. See also Stouffer's method for combining Z scores, which allows weighting.

Usage

FishersMethod(pvalues, ...)

Arguments

Value

• A meta-analytic p-value

References

• Fisher, R.A. (1925). Statistical Methods for Research Workers. Oliver and Boyd (Edinburgh). ISBN 0-05-002170-2.

• Fisher, R. A (1948). "Questions and answers #14". The American Statistician. 2: 30–31. doi:10.2307/2681650.

• Stouffer, S. A. and Suchman, E. A. and DeVinney, L. C. and Star, S. A. and Williams, R. M. Jr. (1949) The American Soldier, Vol. 1 - Adjustment during Army Life. Princeton, Princeton University Press.

See Also

Other Miscellaneous Stats Functions: SE_from_p(), geometric_mean(), harmonic_mean(), oddsratio(), reliability(), umx, umxCov2cor(), umxHetCor(), umxParan(), umxWeightedAIC(), umx_apply(), umx_cor(), umx_means(), umx_r_test(), umx_round(), umx_scale(), umx_var()

Examples

FishersMethod(c(.041, .378))

Twin data: General Family Functioning, divorce, and well-being.

Description

Measures of family functioning, happiness and related variables in twins, and their brothers and sisters. (see details)

Usage

data(GFF)

Format

A data frame with 1000 rows of twin-family data columns.

Details

Several scales in the data are described in van der Aa et al. (2010). General Family Functioning (GFF) refers to adolescents' evaluations general family health vs. pathology. It assesses problem solving, communication, roles within the household, affection, and control. GFF was assessed with a Dutch translation of the General Functioning sub-scale of the McMaster Family Assessment Device (FAD) (Epstein et al., 1983).

Family Conflict (FC) refers to adolescents' evaluations of the amount of openly expressed anger, aggression, and conflict among family members. Conflict sub-scale of the Family Environment Scale (FES) (Moos, 1974)

Quality of life in general (QLg) was assessed with the 10-step Cantril Ladder from best- to worst-possible life (Cantril, 1965).

• zyg_6grp: Six-level zygosity: MZMM, DZMM, MZFF, DZFF, DZMF, DZFM

• zyg_2grp: Two-level zygosity measure: 'MZ', 'DZ'

• divorce: Parental divorce status: 0 = No, 1 = Yes

• sex_T1: Sex of twin 1: 0 = "male", 1 = "female"

• age_T1: Age of twin 1 (years)

• gff_T1: General family functioning for twin 1

• fc_T1: Family conflict sub-scale of the FES

• qol_T1: Quality of life for twin 1

• hap_T1: General happiness for twin 1

• sat_T1: Satisfaction with life for twin 1

• AD_T1: Anxiety and Depression for twin 1

• SOMA_T1: Somatic complaints for twin 1

• SOC_T1: Social problems for twin 1

• THOU_T1: Thought disorder problems for twin 1

• sex_T2: Sex of twin 2

• age_T2: Age of twin 2

• gff_T2: General family functioning for twin 2

• fc_T2: Family conflict sub-scale of the FES

• qol_T2: Quality of life for twin 2

• hap_T2: General happiness for twin 2

• sat_T2: Satisfaction with life for twin 2

• AD_T2: Anxiety and Depression for twin 2

• SOMA_T2: Somatic complaints for twin 2

• SOC_T2: Social problems for twin 2

• THOU_T2: Thought disorder problems for twin 2

• sex_Ta: Sex of sib 1

• age_Ta: Age of sib 1

• gff_Ta: General family functioning for sib 1

• fc_Ta: Family conflict sub-scale of the FES

• qol_Ta: Quality of life for sib 1

• hap_Ta: General happiness for sib 1

• sat_Ta: Satisfaction with life for sib 1

• AD_Ta: Anxiety and Depression for sib 1

• SOMA_Ta: Somatic complaints for sib 1

• SOC_Ta: Social problems for sib 1

• THOU_Ta: Thought disorder problems for sib 1

• sex_Ts: Sex of sib 2

• age_Ts: Age of sib 2

• gff_Ts: General family functioning for sib 2

• fc_Ts: Family conflict sub-scale of the FES

• qol_Ts: Quality of life for sib 2

• hap_Ts: General happiness for sib 2

• sat_Ts: Satisfaction with life for sib 2

• AD_Ts: Anxiety and Depression for sib 2

• SOMA_Ts: Somatic complaints for sib 2

• SOC_Ts: Social problems for sib 2

• THOU_Ts: Thought disorder problems for sib 2

References

van der Aa, N., Boomsma, D. I., Rebollo-Mesa, I., Hudziak, J. J., & Bartels, M. (2010). Moderation of genetic factors by parental divorce in adolescents' evaluations of family functioning and subjective wellbeing. Twin Research and Human Genetics, 13, 143-162. doi:10.1375/twin.13.2.143

See Also

Other datasets: Fischbein_wt, docData, iqdat, umx, us_skinfold_data

Examples

## Not run: 
# Twin 1 variables (end in '_T1')
data(GFF)
umx_names(GFF, "1$") # Just variables ending in 1 (twin 1)
str(GFF) # first few rows


m1 = umxACE(selDVs= "gff", sep = "_T",
	mzData = subset(GFF, zyg_2grp == "MZ"), 
	dzData = subset(GFF, zyg_2grp == "DZ")
)

## End(Not run)

Generic RMSEA function

Description

See RMSEA.MxModel() to access the RMSEA of MxModels

Usage

RMSEA(x, ci.lower, ci.upper, digits)

Arguments

Value

• RMSEA object containing value (and perhaps a CI)

See Also

Other Reporting functions: RMSEA.MxModel(), RMSEA.summary.mxmodel(), extractAIC.MxModel(), loadings(), loadings.MxModel(), residuals.MxModel(), tmx_show(), tmx_show.MxMatrix(), umxCI(), umxCI_boot(), umxConfint(), umxExpCov(), umxExpMeans(), umxFitIndices(), umxRotate()

RMSEA function for MxModels

Description

Return RMSEA and its confidence interval on a model. RMSEA(tmp, silent=TRUE)$RMSEA

Usage

## S3 method for class 'MxModel'
RMSEA(x, ci.lower = 0.025, ci.upper = 0.975, digits = 3)

Arguments

Value

• object containing the RMSEA, lower and upper bounds, and p-close

References

• https://github.com/tbates/umx

See Also

Other Reporting functions: RMSEA(), RMSEA.summary.mxmodel(), extractAIC.MxModel(), loadings(), loadings.MxModel(), residuals.MxModel(), tmx_show(), tmx_show.MxMatrix(), umxCI(), umxCI_boot(), umxConfint(), umxExpCov(), umxExpMeans(), umxFitIndices(), umxRotate()

Examples

## Not run: 
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)
)
RMSEA(m1)

x = RMSEA(m1)
x$RMSEA # 0.0309761

# Raw: needs to be run by umx to get RMSEA
m2 = umxRAM("One Factor", data = demoOneFactor,
	umxPath("G", to = manifests),
	umxPath(v.m. = manifests),
	umxPath(v1m0 = "G")
)
RMSEA(m2)

## End(Not run)

RMSEA function for MxModel summary

Description

Extract the RMSEA and confidence interval from a model summary and returns it as an RMSEA object. To report just the RMSEA, you can use RMSEA(model)$RMSEA

Usage

## S3 method for class 'summary.mxmodel'
RMSEA(x, ci.lower = 0.025, ci.upper = 0.975, digits = 3)

Arguments

Value

• object containing the RMSEA and lower and upper bounds

References

• https://github.com/simsem/semTools/wiki/Functions, https://github.com/tbates/umx

See Also

Other Reporting functions: RMSEA(), RMSEA.MxModel(), extractAIC.MxModel(), loadings(), loadings.MxModel(), residuals.MxModel(), tmx_show(), tmx_show.MxMatrix(), umxCI(), umxCI_boot(), umxConfint(), umxExpCov(), umxExpMeans(), umxFitIndices(), umxRotate()

Examples

## Not run: 
require(umx)
data(demoOneFactor)
manifests = names(demoOneFactor)

m1 = umxRAM("One Factor", data = demoOneFactor[1:100,], type = "cov",
	umxPath("G", to = manifests),
	umxPath(var = manifests),
	umxPath(var = "G", fixedAt = 1.0)
)
tmp = summary(m1)
RMSEA(tmp)

## End(Not run)

Compute an SE from a beta and p value

Description

SE_from_p takes beta and p, and returns an SE.

Usage

SE_from_p(beta = NULL, p = NULL, SE = NULL, lower = NULL, upper = NULL)

Arguments

Value

• Standard error

See Also

• umxAPA()

Other Miscellaneous Stats Functions: FishersMethod(), geometric_mean(), harmonic_mean(), oddsratio(), reliability(), umx, umxCov2cor(), umxHetCor(), umxParan(), umxWeightedAIC(), umx_apply(), umx_cor(), umx_means(), umx_r_test(), umx_round(), umx_scale(), umx_var()

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)

Print a money object

Description

Print function for "money" objects, e.g. fin_interest().

Usage

bucks(
  x,
  symbol = "$",
  big.mark = ",",
  decimal.mark = ".",
  trim = TRUE,
  largest_with_cents = 1e+05,
  negative_parens = c("hyphen", "minus", "parens"),
  ...
)

Arguments

Value

• invisible

See Also

• fin_percent(), fin_interest(), scales::dollar()

Examples

bucks(100 * 1.05^32)
fin_interest(deposits = 20e3, interest = 0.07, yrs = 20)

Convert Degrees to Degrees

Description

A helper to convert degrees (360 in a circle) to Rad (2\pi in a circle).

note: R's trig functions, e.g. sin() use Radians for input!

The formula is radians = deg x 180/\pi.

• 180 Degrees is equal to \pi radians.

• 1 Rad = 180/\pi degrees = ~ 57.296 degrees.

Usage

deg2rad(deg)

Arguments

Value

• value in radians

References

https://en.wikipedia.org/wiki/Radian

See Also

• rad2deg(), sin()

Other Miscellaneous Functions: fin_JustifiedPE(), fin_NI(), fin_interest(), fin_percent(), fin_ticker(), fin_valuation(), rad2deg(), umxBrownie()

Examples

deg2rad(180) == pi # TRUE!

dl_from_dropbox

Description

Download a file from Dropbox, given either the url, or the name and key

Usage

dl_from_dropbox(x, key = NULL)

Arguments

Details

Improvements would include error handling...

Value

None

References

- https://thebiobucket.blogspot.kr/2013/04/download-files-from-dropbox.html

See Also

Other File Functions: umx, umx_file_load_pseudo(), umx_make_sql_from_excel(), umx_move_file(), umx_open(), umx_rename_file(), umx_write_to_clipboard()

Examples

## Not run: 
dl_from_dropbox("https://dl.dropboxusercontent.com/s/7kauod48r9cfhwc/tinytwinData.rda")
dl_from_dropbox("tinytwinData.rda", key = "7kauod48r9cfhwc")

## End(Not run)

Twin data for Direction of causation modelling

Description

A dataset containing indicators for two traits varA and varB, each measured in MZ and DZ twins.

Usage

data(docData)

Format

A data frame 6 manifests for each of two twins in 1400 families of MZ and DZ twins

Details

It is designed to show off umxDoC() testing the hypothesis varA causes varB, varB causes varA, both cause each other.

• zygosity "MZFF", "DZFF", "MZMM", or "DZMM"

• varA1_T1 Twin one's manifest 1 for varA

• varA2_T1 Twin one's manifest 2 for varA

• varA3_T1 Twin one's manifest 3 for varA

• varB1_T1 Twin one's manifest 1 for varB

• varB2_T1 Twin one's manifest 2 for varB

• varB3_T1 Twin one's manifest 3 for varB

• varA1_T2 Twin two's manifest 1 for varA

• varA2_T2 Twin two's manifest 2 for varA

• varA3_T2 Twin two's manifest 3 for varA

• varB1_T2 Twin two's manifest 1 for varB

• varB2_T2 Twin two's manifest 2 for varB

• varB3_T2 Twin two's manifest 3 for varB

References

• N.A. Gillespie and N.G. Martin (2005). Direction of Causation Models. In Encyclopedia of Statistics in Behavioral Science, 1, 496–499. Eds. Brian S. Everitt & David C. Howell

See Also

• umxDoC(), plot.MxModelDoC(), umxSummary.MxModelDoC(), umxModify()

Other datasets: Fischbein_wt, GFF, iqdat, umx, us_skinfold_data

Examples

data(docData)
str(docData)
mzData = subset(docData, zygosity %in% c("MZFF", "MZMM"))
dzData = subset(docData, zygosity %in% c("DZFF", "DZMM"))
par(mfrow = c(1, 2))  # 1 rows and 3 columns
plot(varA1_T2 ~varA1_T1, ylim = c(-4, 4), data = mzData, main="MZ")
tmp = round(cor.test(~varA1_T1 + varA1_T2, data = mzData)$estimate, 2)
text(x=-4, y=3, labels = paste0("r = ", tmp))
plot(varA1_T2 ~varA1_T1, ylim = c(-4, 4), data = dzData, main="DZ")
tmp = round(cor.test(~varA1_T1 + varA1_T2, data = dzData)$estimate, 2)
text(x=-4, y=3, labels = paste0("r = ", tmp))
par(mfrow = c(1, 1))  # back to as it was

Extract AIC from MxModel

Description

Returns the AIC for an OpenMx model. Original Author: Brandmaier

Usage

## S3 method for class 'MxModel'
extractAIC(fit, scale, k, ...)

Arguments

Value

• AIC value

See Also

• AIC(), umxCompare(), logLik()

Other Reporting functions: RMSEA(), RMSEA.MxModel(), RMSEA.summary.mxmodel(), loadings(), loadings.MxModel(), residuals.MxModel(), tmx_show(), tmx_show.MxMatrix(), umxCI(), umxCI_boot(), umxConfint(), umxExpCov(), umxExpMeans(), umxFitIndices(), umxRotate()

Examples

## Not run: 
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)
)
extractAIC(m1)
# -2.615998
AIC(m1)

## End(Not run)

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.

Usage

fin_JustifiedPE(
  Dividend = 0.02,
  EPS = 1,
  growthRate = 0.08,
  discountRate = 0.12,
  basePE = 20,
  yrs = 10
)

Arguments

Value

• A PE that is justified for this stock.

See Also

• fin_interest(), fin_percent(), fin_NI()

Other Miscellaneous Functions: deg2rad(), fin_NI(), fin_interest(), fin_percent(), fin_ticker(), fin_valuation(), rad2deg(), umxBrownie()

Examples

# fin_JustifiedPE(Dividend= .8, EPS = 2, growthRate = .06, discountRate = .1)

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.

Usage

fin_NI(annualEarnings, symbol = "£")

Arguments

Value

• NI

References

• https://www.telegraph.co.uk/tax/tax-hacks/politicians-running-scared-long-overdue-national-insurance-overhaul/

See Also

• fin_interest(), fin_percent(), fin_valuation()

Other Miscellaneous Functions: deg2rad(), fin_JustifiedPE(), fin_interest(), fin_percent(), fin_ticker(), fin_valuation(), rad2deg(), umxBrownie()

Examples

fin_NI(42e3)
fin_NI(142000)

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+.53 = 10.5 up to 9+210 = 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)

Usage

fin_interest(
  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")
)

Arguments

Value

• Value of balance after yrs of investment.

References

• https://en.wikipedia.org/wiki/Compound_interest

See Also

• umx_set_dollar_symbol(), fin_percent(), fin_NI(), fin_valuation()

Other Miscellaneous Functions: deg2rad(), fin_JustifiedPE(), fin_NI(), fin_percent(), fin_ticker(), fin_valuation(), rad2deg(), umxBrownie()

Examples

## Not run: 
# 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")

## End(Not run)

# 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)

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%)

Usage

fin_percent(
  percent,
  value = 100,
  symbol = "$",
  digits = 2,
  plot = TRUE,
  logY = TRUE
)

Arguments

Value

• new value and change required to return to baseline.

See Also

• fin_interest()

Other Miscellaneous Functions: deg2rad(), fin_JustifiedPE(), fin_NI(), fin_interest(), fin_ticker(), fin_valuation(), rad2deg(), umxBrownie()

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)

Open a ticker in yahoo finance.

Description

Open a stock ticker, currently in yahoo finance

Usage

fin_ticker(ticker = "INTC")

Arguments

Value

• Open a ticker in a finance site online

See Also

• fin_interest(), fin_percent(), fin_NI()

Other Miscellaneous Functions: deg2rad(), fin_JustifiedPE(), fin_NI(), fin_interest(), fin_percent(), fin_valuation(), rad2deg(), umxBrownie()

Examples

# Open $INTC in yahoo finance.
## Not run: 
fin_ticker("INTC")

## End(Not run)

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.

Usage

fin_valuation(
  revenue = 6e+06 * 30000,
  opmargin = 0.08,
  expenses = 0.2,
  PE = 30,
  symbol = "$",
  use = c("B", "M")
)

Arguments

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

Value

• value

See Also

• fin_interest(), fin_NI(), fin_percent()

Other Miscellaneous Functions: deg2rad(), fin_JustifiedPE(), fin_NI(), fin_interest(), fin_percent(), fin_ticker(), rad2deg(), umxBrownie()

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

Geometric Mean

Description

The Geometric mean is the nth-root of the product of n input values. Common uses include computing economic utility. For example, the geometric mean utility of c(1, 2, 10) is

(1 * 2 * 10)^\frac{1}{3}

= 2.7 not 4.3 (the arithmetic mean of utility).

Usage

geometric_mean(x, na.rm = c(TRUE, FALSE))

Arguments

Value

• Geometric mean of x

References

• https://en.wikipedia.org/wiki/Geometric_mean

See Also

• harmonic_mean(), mean()

Other Miscellaneous Stats Functions: FishersMethod(), SE_from_p(), harmonic_mean(), oddsratio(), reliability(), umx, umxCov2cor(), umxHetCor(), umxParan(), umxWeightedAIC(), umx_apply(), umx_cor(), umx_means(), umx_r_test(), umx_round(), umx_scale(), umx_var()

Examples

geometric_mean(c(50, 100))

# For a given sum, geometric mean is maximised when all values are equal:
geometric_mean(c(75,75))

v = c(1, 149); c(sum(v), geometric_mean(v), mean(v), median(v))
# 150.00000  12.20656  75.00000  75.00000

# Underlying logic
sqrt(50 * 100)

# Alternate form using logs
exp(mean(log(c(50 *100))))

# Reciprocal duality
1/geometric_mean(c(100, 50))
geometric_mean(c(1/100, 1/50))

Add a fit statistic to a ggplot

Description

Add a fit statistic to a ggplot

Usage

ggAddR(model, effect = NA, xloc = 8, yloc = 10)

Arguments

Value

• plot

See Also

• umxPlot(), umxPlotFun()

Other Plotting functions: plot.MxLISRELModel(), plot.MxModel(), plot.MxModelTwinMaker(), umx, umxPlot(), umxPlotACE(), umxPlotACEcov(), umxPlotACEv(), umxPlotCP(), umxPlotDoC(), umxPlotFun(), umxPlotGxE(), umxPlotGxEbiv(), umxPlotIP(), umxPlotSexLim(), umxPlotSimplex()

Examples

## Not run: 
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

## End(Not run)

Harmonic Mean

Description

The harmonic mean is the reciprocal of the arithmetic mean of the reciprocals of the input values. Common uses include computing the mean of ratios, for instance the average P/E ratio in a portfolio. Also it is the correct mean for averaging speeds weighted for distance.

Usage

harmonic_mean(x, weights = NULL, na.rm = c(TRUE, FALSE))

Arguments

Value

• Harmonic mean of x

References

• https://en.wikipedia.org/wiki/Harmonic_mean

See Also

• geometric_mean(), aggregate()

Other Miscellaneous Stats Functions: FishersMethod(), SE_from_p(), geometric_mean(), oddsratio(), reliability(), umx, umxCov2cor(), umxHetCor(), umxParan(), umxWeightedAIC(), umx_apply(), umx_cor(), umx_means(), umx_r_test(), umx_round(), umx_scale(), umx_var()

Examples

# Harmonic means are suitable for ratios
tmp = c(33/1, 23/1)
harmonic_mean(tmp)

geometric_mean(tmp)
mean(tmp)

# Example with weights
harmonic_mean(c(33/1, 23/1), weights= c(.2, .8))
# If Jack travels outbound at 1 mph, and returns at 10 miles an hour, what is his average speed?
harmonic_mean(c(1,10)) # 1.81 mph

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.

Usage

install.OpenMx(
  loc = c("NPSOL", "travis", "CRAN", "open travis build page", "UVa"),
  url = NULL,
  lib,
  repos = getOption("repos")
)

Arguments

Value

None

References

• https://github.com/tbates/umx, https://tbates.github.io

See Also

umxVersion()

Other Miscellaneous Utility Functions: libs(), qm(), umx, umxLav2RAM(), umxModelNames(), umxRAM2Lav(), umxVersion(), umx_array_shift(), umx_find_object(), umx_lower.tri(), umx_msg(), umx_open_CRAN_page(), umx_pad(), umx_print()

Examples

## Not run: 
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

## End(Not run)

Twin data: IQ measured longitudinally across 4 ages.

Description

Measures of IQ across four ages in 261 pairs of identical twins and 301 pairs of fraternal (DZ) twins. (see details). It is used as data for the [umxSimplex()] examples.

Usage

data(iqdat)

Format

A data frame with 562 rows (twin families). Nine measures on each twin.

Details

• zygosity Zygosity (MZ or DZ)

• IQ_age1_T1 T1 IQ measured at age 1

• IQ_age2_T1 T1 IQ measured at age 2

• IQ_age3_T1 T1 IQ measured at age 3

• IQ_age4_T1 T1 IQ measured at age 4

• IQ_age1_T2 T2 IQ measured at age 1

• IQ_age2_T2 T2 IQ measured at age 2

• IQ_age3_T2 T2 IQ measured at age 3

• IQ_age4_T2 T2 IQ measured at age 4

References

Boomsma, D. I., Martin, N. G., & Molenaar, P. C. (1989). Factor and simplex models for repeated measures: application to two psychomotor measures of alcohol sensitivity in twins. *Behavior Genetics*, **19**, 79-96. Retrieved from <https://www.ncbi.nlm.nih.gov/pubmed/2712815>

See Also

[umxSimplex()]

Other datasets: Fischbein_wt, GFF, docData, umx, us_skinfold_data

Examples

## Not run: 
data(iqdat)
str(iqdat)
par(mfrow = c(1, 3))  # 1 rows and 3 columns
plot(IQ_age4_T1 ~ IQ_age4_T2, ylim = c(50, 150), data = subset(iqdat, zygosity == "MZ"))
plot(IQ_age4_T1 ~ IQ_age4_T2, ylim = c(50, 150), data = subset(iqdat, zygosity == "DZ"))
plot(IQ_age1_T1 ~ IQ_age4_T2, data = subset(iqdat, zygosity == "MZ"))
par(mfrow = c(1, 1))  # back to as it was

## End(Not run)

load libraries

Description

libs allows loading multiple libraries in one call

Usage

libs(..., force.update = FALSE)

Arguments

Value

• nothing.

See Also

• library(), install.packages(), remove.packages()

Other Miscellaneous Utility Functions: install.OpenMx(), qm(), umx, umxLav2RAM(), umxModelNames(), umxRAM2Lav(), umxVersion(), umx_array_shift(), umx_find_object(), umx_lower.tri(), umx_msg(), umx_open_CRAN_page(), umx_pad(), umx_print()

Examples

## Not run: 
libs("umx", "OpenMx", "car")
libs("umx", c("OpenMx", "car"))
remove.packages()

## End(Not run)

loadings Generic loadings function to extract factor loadings from exploratory or confirmatory factor analyses.

Description

See loadings.MxModel to access the loadings of OpenMx EFA models.

Usage

loadings(x, ...)

Arguments

Details

Base loadings handles factanal() objects.

Value

• matrix of loadings

References

• https://tbates.github.io, https://github.com/tbates/umx

See Also

Other Reporting functions: RMSEA(), RMSEA.MxModel(), RMSEA.summary.mxmodel(), extractAIC.MxModel(), loadings.MxModel(), residuals.MxModel(), tmx_show(), tmx_show.MxMatrix(), umxCI(), umxCI_boot(), umxConfint(), umxExpCov(), umxExpMeans(), umxFitIndices(), umxRotate()

Extract factor loadings from an EFA (factor analysis).

Description

loadings extracts the factor loadings from an EFA (factor analysis) model. It behaves equivalently to stats::loadings, returning the loadings from an EFA (factor analysis). However it does not store the rotation matrix.

Usage

## S3 method for class 'MxModel'
loadings(x, ...)

Arguments

Value

• loadings matrix

References

• https://github.com/tbates/umx, https://tbates.github.io

See Also

• factanal(), loadings()

Other Reporting functions: RMSEA(), RMSEA.MxModel(), RMSEA.summary.mxmodel(), extractAIC.MxModel(), loadings(), residuals.MxModel(), tmx_show(), tmx_show.MxMatrix(), umxCI(), umxCI_boot(), umxConfint(), umxExpCov(), umxExpMeans(), umxFitIndices(), umxRotate()

Examples

## Not run: 
myVars = c("mpg", "disp", "hp", "wt", "qsec")
m1 = umxEFA(name = "test", factors = 2, data = mtcars[, myVars])
loadings(m1)

## End(Not run)

Succinctly select complete rows from a dataframe

Description

Succinctly select complete rows from a dataframe.

Usage

noNAs(df, rows = NULL, cols = NULL, drop = TRUE)

Arguments

Value

• Complete rows and (optionally) selected columns

See Also

Other Data Functions: prolific_anonymize(), prolific_check_ID(), prolific_read_demog(), umx, umxFactor(), umxHetCor(), umx_as_numeric(), umx_cont_2_quantiles(), umx_lower2full(), umx_make_MR_data(), umx_make_TwinData(), umx_make_fake_data(), umx_make_raw_from_cov(), umx_merge_randomized_columns(), umx_polychoric(), umx_polypairwise(), umx_polytriowise(), umx_read_lower(), umx_rename(), umx_reorder(), umx_score_scale(), umx_select_valid(), umx_stack(), umx_strings2numeric()

Examples

tmp = mtcars
tmp[2,1] = NA
noNAs(tmp, cols="mpg")
noNAs(tmp, cols="mpg", drop = FALSE)
noNAs(tmp) # no Mazda RX4 Wag

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.

Usage

oddsratio(grp1 = c(n = 3, N = 10), grp2 = c(n = 1, N = 10), alpha = 0.05)

Arguments

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]

Value

• List of odds in group 1 and group2, and the resulting OR and CI

References

• https://stats.oarc.ucla.edu/r/dae/logit-regression/, https://tbates.github.io

See Also

• umx_r_test()

Other Miscellaneous Stats Functions: FishersMethod(), SE_from_p(), geometric_mean(), harmonic_mean(), reliability(), umx, umxCov2cor(), umxHetCor(), umxParan(), umxWeightedAIC(), umx_apply(), umx_cor(), umx_means(), umx_r_test(), umx_round(), umx_scale(), umx_var()

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)

Create and display a graphical path diagram for a LISREL model.

Description

plot.MxLISRELModel produces SEM diagrams using [DiagrammeR::DiagrammeR() to create the image.

Usage

## S3 method for class 'MxLISRELModel'
plot(
  x = NA,
  std = FALSE,
  fixed = TRUE,
  means = TRUE,
  digits = 2,
  file = "name",
  labels = c("none", "labels", "both"),
  resid = c("circle", "line", "none"),
  strip_zero = TRUE,
  ...
)

Arguments

Details

Note: By default, plots open in your browser (or plot pane if using RStudio).

Opening in an external editor/app

The underlying format is graphviz. If you use umx_set_plot_format("graphviz"), figures will open in a graphviz helper app (if installed). If you use graphviz, we try and use that app, but YOU HAVE TO INSTALL IT!

On MacOS, you may need to associate the ‘.gv’ extension with your graphviz app. Find the .gv file made by plot, get info (cmd-I), then choose “open with”, select graphviz.app (or OmniGraffle professional), then set “change all”.

The commercial application “OmniGraffle” is great for editing these images.

References

• https://github.com/tbates/umx, https://en.wikipedia.org/wiki/DOT_(graph_description_language)

See Also

• umx_set_plot_format(), umx_set_auto_plot(), umx_set_plot_format(), plot.MxModel(), umxPlotACE(), umxPlotCP(), umxPlotIP(), umxPlotGxE()

Other Plotting functions: ggAddR(), plot.MxModel(), plot.MxModelTwinMaker(), umx, umxPlot(), umxPlotACE(), umxPlotACEcov(), umxPlotACEv(), umxPlotCP(), umxPlotDoC(), umxPlotFun(), umxPlotGxE(), umxPlotGxEbiv(), umxPlotIP(), umxPlotSexLim(), umxPlotSimplex()

Examples

# plot()
# TODO get LISREL example model
# Figure out how to map its matrices to plot. Don't do without establishing demand.

Create and display a graphical path diagram for a model.

Description

plot() produces SEM diagrams in graphviz format, and relies on DiagrammeR::DiagrammeR() to create the image.

Usage

## S3 method for class 'MxModel'
plot(
  x = NA,
  std = FALSE,
  fixed = TRUE,
  means = TRUE,
  digits = 2,
  file = "name",
  labels = c("none", "labels", "both"),
  resid = c("circle", "line", "none"),
  strip_zero = FALSE,
  splines = c("TRUE", "FALSE", "compound", "ortho", "polyline"),
  min = NULL,
  same = NULL,
  max = NULL,
  ...
)

Arguments

Details

Note: DiagrammeR::DiagrammeR() is supported out of the box. By default, plots open in your browser. Other options include pdf SVG etc.

If you use umx_set_plot_format("graphviz"), graphs will open in a graphviz helper app (if installed).

The commercial application “OmniGraffle” is great for editing these images. On unix and windows, plot() will create a pdf and open it in your default pdf reader.

If you use graphviz, we try and use that app, but YOU HAVE TO INSTALL IT!

MacOS note: On Mac, we will try and open an app: you may need to associate the ‘.gv’ extension with the graphviz app. Find the .gv file made by plot, get info (cmd-I), then choose “open with”, select graphviz.app (or OmniGraffle professional), then set “change all”.

References

• https://github.com/tbates/umx, https://en.wikipedia.org/wiki/DOT_(graph_description_language)

See Also

• umx_set_plot_format(), plot.MxModel(), umxPlotACE(), umxPlotCP(), umxPlotIP(), umxPlotGxE()

Other Plotting functions: ggAddR(), plot.MxLISRELModel(), plot.MxModelTwinMaker(), umx, umxPlot(), umxPlotACE(), umxPlotACEcov(), umxPlotACEv(), umxPlotCP(), umxPlotDoC(), umxPlotFun(), umxPlotGxE(), umxPlotGxEbiv(), umxPlotIP(), umxPlotSexLim(), umxPlotSimplex()

Examples

## Not run: 
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)
)
plot(m1)
plot(m1, std = TRUE, resid = "line", digits = 3, strip_zero = FALSE)

# ============================================================
# = With a growth model, demonstrate splines= false to force =
# = straight lines, and move "rank" of intercept object      =
# ============================================================

m1 = umxRAM("grow", data = myGrowthMixtureData,
	umxPath(var = manifests, free = TRUE), 
	umxPath(means = manifests, fixedAt = 0), 
	umxPath(v.m. = c("int","slope")),
	umxPath("int", with = "slope"),
	umxPath("int", to = manifests, fixedAt = 1), 
	umxPath("slope", to = manifests, arrows = 1, fixedAt = c(0,1,2,3,4))
)

plot(m1, means=FALSE, strip=TRUE, splines="FALSE", max="int")

## End(Not run) # end dontrun

Create and display a graphical path diagram for a path-based twin model.

Description

Assumes the model has a group called "MZ" inside.

Usage

## S3 method for class 'MxModelTwinMaker'
plot(
  x = NA,
  std = FALSE,
  fixed = TRUE,
  means = TRUE,
  oneTwin = TRUE,
  sep = "_T",
  digits = 2,
  file = "name",
  labels = c("none", "labels", "both"),
  resid = c("circle", "line", "none"),
  strip_zero = FALSE,
  splines = TRUE,
  min = NULL,
  same = NULL,
  max = NULL,
  ...
)

Arguments

Details

If you use umx_set_plot_format("graphviz"), they will open in a graphviz helper app (if installed). The commercial application “OmniGraffle” is great for editing these images. On unix and windows, plot() will create a pdf and open it in your default pdf reader.

See Also

• umx_set_plot_format(), plot.MxModel(), umxPlotACE(), umxPlotCP(), umxPlotIP(), umxPlotGxE()

Other Plotting functions: ggAddR(), plot.MxLISRELModel(), plot.MxModel(), umx, umxPlot(), umxPlotACE(), umxPlotACEcov(), umxPlotACEv(), umxPlotCP(), umxPlotDoC(), umxPlotFun(), umxPlotGxE(), umxPlotGxEbiv(), umxPlotIP(), umxPlotSexLim(), umxPlotSimplex()

Examples

## Not run: 
require(umx)
# 
# =====================
# = Make an ACE model =
# =====================
# 1. Clean data: Add separator and scale
data(twinData)
tmp = umx_make_twin_data_nice(data=twinData, sep="", zygosity="zygosity", numbering=1:2)
tmp = umx_scale_wide_twin_data(varsToScale= c("wt", "ht"), sep= "_T", data= tmp)
mzData = subset(tmp, zygosity %in%  c("MZFF", "MZMM"))
dzData = subset(tmp, zygosity %in%  c("DZFF", "DZMM"))

# 2. Define paths: You only need the paths for one person:
paths = c(
umxPath(v1m0 = c("a1", 'c1', "e1")),
umxPath(means = c("wt")),
umxPath(c("a1", 'c1', "e1"), to = "wt", values=.2)
)
m1 = umxTwinMaker("test", paths, mzData = mzData, dzData= dzData)
plot(m1, std= TRUE, means= FALSE)
plot(m1, means=FALSE, std=TRUE, strip=TRUE, splines="FALSE", max="intercept")

## End(Not run) # end dontrun

# =================
# = An ACEv model =
# =================
# Not complete

paths = c(
umxPath(v1m0 = c("A1", 'C1', "E1")),
umxPath(v1m0 = c("A2", 'C2', "E2")),
umxPath(v.m0 = c("l1", 'l2')),
umxPath(v.m. = c("wt", "ht")),
umxPath(c("A1", 'C1', "E1"), to = "l1", values= .2),
umxPath(c("A2", 'C2', "E2"), to = "l2", values= .2),
umxPath(c("l1", 'l2'), to = c("wt", "ht"), values= .2)
)

Plot a percent change graph

Description

Plot method for "percent" objects: e.g. fin_percent().

Usage

## S3 method for class 'percent'
plot(x, ...)

Arguments

Value

• invisible

See Also

• fin_percent()

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)

Test the power of an ACE model to detect paths of interest.

Description

power.ACE.test simulates a univariate ACE model. It computes power to detect dropping one or more paths (a, c, or a after dropping c), specified in ⁠drop=⁠.

The interface and functionality of this service are experimental and subject to change.

Usage

power.ACE.test(
  AA = 0.5,
  CC = 0,
  EE = NULL,
  DD = NULL,
  update = c("a", "c", "a_after_dropping_c", "d"),
  value = 0,
  n = NULL,
  MZ_DZ_ratio = 1,
  sig.level = 0.05,
  power = 0.8,
  method = c("ncp", "empirical"),
  search = FALSE,
  tryHard = c("yes", "no", "ordinal", "search"),
  digits = 2,
  optimizer = NULL,
  nSim = 4000
)

Arguments

Details

Statistical power is the proportion of studies that, over the long run, one should expect to yield a statistically significant result given certain study characteristics such as sample size (N), the expected effect size (\beta), and the criterion for statistical significance (\alpha).

(with nMZpairs= 2000 and MZ_DZ_ratio*nMZpairs DZ twins.

A typical target for power is 80%. Much as the accepted critical p-value is .05, this has emerged as a trade off, in this case of resources required for more powerful studies against the cost of missing a true effect. People interested in truth discourage running studies with low power: A study with 20 percent power will fail to detect real effects 80% of the time. But even with zero power, the Type-I error rate remains a nominal 5% (and with any researcher degrees of freedom, perhaps much more than that). Low powered research, then, fails to detect true effects, and generates support for random false theories about as often. This sounds silly, but empirical rates are often as low as 20% (Button, et al., 2013).

Illustration of \alpha, \beta, and power (1-\beta):

Value

OpenMx::mxPower() object

References

• Visscher, P.M., Gordon, S., Neale, M.C. (2008). Power of the classical twin design revisited: II detection of common environmental variance. Twin Res Hum Genet, 11: 48-54. doi:10.1375/twin.11.1.48.

• Button, K. S., Ioannidis, J. P., Mokrysz, C., Nosek, B. A., Flint, J., Robinson, E. S., and Munafo, M. R. (2013). Power failure: why small sample size undermines the reliability of neuroscience. Nature Reviews Neuroscience, 14, 365-376. doi:10.1038/nrn3475

See Also

• umxPower(), OpenMx::mxPower(), umxACE()

Other Twin Modeling Functions: umx, umxACE(), umxACEcov(), umxACEv(), umxCP(), umxDiffMZ(), umxDiscTwin(), umxDoC(), umxDoCp(), umxGxE(), umxGxE_window(), umxGxEbiv(), umxIP(), umxMRDoC(), umxReduce(), umxReduceACE(), umxReduceGxE(), umxRotate.MxModelCP(), umxSexLim(), umxSimplex(), umxSummarizeTwinData(), umxSummaryACE(), umxSummaryACEv(), umxSummaryDoC(), umxSummaryGxEbiv(), umxSummarySexLim(), umxSummarySimplex(), umxTwinMaker()

Examples

# =====================================================
# = N for .8 power to detect a^2 = .5 equal MZ and DZ =
# =====================================================
power.ACE.test(AA = .5, CC = 0, update = "a")
# Suggests n = 84 MZ and 94 DZ pairs.

## Not run: 
# ================================
# = Show power across range of N =
# ================================
power.ACE.test(AA= .5, CC= 0, update = "a", search = TRUE)

# Salutary note: You need well fitting models with correct betas in the data
# for power to be valid.
# tryHard helps ensure this, as does the default nSim= 4000 pair data.
# Power is important to get right, so I recommend using tryHard = "yes" (the default)

# =====================
# = Power to detect C =
# =====================

# 102 of each of MZ and DZ pairs for 80% power (default).
power.ACE.test(AA= .5, CC= .3, update = "c")

# ==========================================
# = Set 'a' to a fixed, but non-zero value =
# ==========================================

power.ACE.test(update= "a", value= sqrt(.2), AA= .5, CC= 0)

# ========================================
# = Drop More than one parameter (A & C) =
# ========================================
# E vs AE: the hypothesis that twins show no familial similarity.
power.ACE.test(update = "a_after_dropping_c", AA= .5, CC= .3)

# ===================================================
# = More power to detect A > 0 when more C present  =
# ===================================================

power.ACE.test(update = "a", AA= .5, CC= .0)
power.ACE.test(update = "a", AA= .5, CC= .3)

# ====================================================
# = More power to detect C > 0 when more A present?  =
# ====================================================

power.ACE.test(update = "c", AA= .0, CC= .5)
power.ACE.test(update = "c", AA= .3, CC= .5)


# ===================================
# = Power with more DZs or more MZs =
# ===================================

# Power about the same: total pairs with 2 MZs per DZ
power.ACE.test(MZ_DZ_ratio= 2/1, update= "a", AA= .3, CC= 0, method="ncp", tryHard="yes")
power.ACE.test(MZ_DZ_ratio= 1/2, update= "a", AA= .3, CC= 0, method="ncp", tryHard="yes")
power.ACE.test(update= "a", AA= .3, CC= 0, method="ncp", tryHard="yes")


# =====================================
# = Compare ncp and empirical methods =
# =====================================

power.ACE.test(update= "a", AA= .5, CC= 0, method = "ncp")
# method = "ncp": For 80% power, you need 166 MZ and 166 DZ pairs
power.ACE.test(update= "a", AA= .5, CC= 0, method= "empirical")
# method= "empirical": For 80% power, you need 154 MZ and 154 DZ pairs

# ====================
# = Show off options =
# ====================
# 1. tryHard

power.ACE.test(update = "a", AA= .5, CC= 0, tryHard= "no")

# 2. toggle optimizer

power.ACE.test(update= "a", AA= .5, CC= 0, optimizer= "SLSQP")

# 3. You can raise or lower the number of pairs used in the true model
#    by varying nSim (twin pairs in the simulated data).

power.ACE.test(update = "a", AA= .5, CC= 0, nSim= 20)


## End(Not run)

Print a RMSEA object

Description

Print method for "RMSEA" objects: e.g. RMSEA().

Usage

## S3 method for class 'RMSEA'
print(x, ...)

Arguments

Value

• invisible

See Also

• RMSEA(), print()

Examples

## Not run: 
data(demoOneFactor)
manifests = names(demoOneFactor)

m1 = umxRAM("One Factor", data = demoOneFactor, type= "cov",
	umxPath("G", to = manifests),
	umxPath(var = manifests),
	umxPath(var = "G", fixedAt = 1.0)
)
tmp = summary(m1)
RMSEA(tmp)

## End(Not run)

Print a scale "oddsratio" object

Description

Print method for the oddsratio() function.

Usage

## S3 method for class 'oddsratio'
print(x, digits = 3, ...)

Arguments

Value

• invisible oddsratio object (x).

See Also

• print(), oddsratio(),

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 a percent object

Description

Print method for "percent" objects: e.g. fin_percent().

Usage

## S3 method for class 'percent'
print(x, ...)

Arguments

Value

• invisible

See Also

• fin_percent()

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 a scale "reliability" object

Description

Print method for the reliability() function.

Usage

## S3 method for class 'reliability'
print(x, digits = 4, ...)

Arguments

Value

• invisible reliability object (x)

See Also

• print(), reliability(),

Examples

# treat vehicle aspects as items of a test
data(mtcars)
reliability(cov(mtcars))

Clean up a prolific file for sharing by removing anonymity-compromising columns.

Description

prolific.ac IDs and other columns like IP and lat/long might compromise subject anonymity when shared. prolific_anonymize replaces PIDs with a simple numeric sequence, preserving repeated measures in long data, and removing other columns. You can delete additional columns by adding them to extraColumns. It is ideal for use when sharing data to https://researchbox.org which enforces anonymization.

Usage

prolific_anonymize(
  df = NULL,
  PID = "PID",
  extraColumns = NA,
  baseOffset = 10000
)

Arguments

Value

• [data.frame]

References

• https://github.com/tbates/umx, https://tbates.github.io

See Also

• prolific_check_ID(), prolific_read_demog(), umx_merge_randomized_columns()

Other Data Functions: noNAs(), prolific_check_ID(), prolific_read_demog(), umx, umxFactor(), umxHetCor(), umx_as_numeric(), umx_cont_2_quantiles(), umx_lower2full(), umx_make_MR_data(), umx_make_TwinData(), umx_make_fake_data(), umx_make_raw_from_cov(), umx_merge_randomized_columns(), umx_polychoric(), umx_polypairwise(), umx_polytriowise(), umx_read_lower(), umx_rename(), umx_reorder(), umx_score_scale(), umx_select_valid(), umx_stack(), umx_strings2numeric()

Examples

## Not run: 
tmp = prolific_anonymize(df, PID = "PID")

## End(Not run)

Return PIDs in df

Description

Participants may time-out on Prolific, but still complete on Qualtrics. This identifies them.

Usage

prolific_check_ID(IDs, df, IDcol = "PROLIFIC_PID")

Arguments

Value

- list of IDs in the dataframe

See Also

- [prolific_read_demog()], [prolific_anonymize()], [umx_merge_randomized_columns()] # [prolific_check_ID()]

Other Data Functions: noNAs(), prolific_anonymize(), prolific_read_demog(), umx, umxFactor(), umxHetCor(), umx_as_numeric(), umx_cont_2_quantiles(), umx_lower2full(), umx_make_MR_data(), umx_make_TwinData(), umx_make_fake_data(), umx_make_raw_from_cov(), umx_merge_randomized_columns(), umx_polychoric(), umx_polypairwise(), umx_polytriowise(), umx_read_lower(), umx_rename(), umx_reorder(), umx_score_scale(), umx_select_valid(), umx_stack(), umx_strings2numeric()

Examples

# IDs = c("59d0ec2446447f00011edb063","5a08c9a7f2e3460001edb063f0254")
# prolific_check_ID(IDs, df)

Read and optionally merge demographics file from prolific academic

Description

prolific academic provides a demographics file. This reads it and merges with your data using PID and participant_id

Usage

prolific_read_demog(
  file,
  base = "",
  df = NULL,
  by.df = "PROLIFIC_PID",
  by.demog = "Participant.id",
  age = "age",
  sex = "Gender",
  vars = NULL,
  all.df = TRUE,
  all.demog = FALSE,
  verbose = FALSE
)

Arguments

Value

• [data.frame]

References

• https://github.com/tbates/umx, https://tbates.github.io

See Also

• prolific_check_ID(), prolific_anonymize(), umx_merge_randomized_columns()

Other Data Functions: noNAs(), prolific_anonymize(), prolific_check_ID(), umx, umxFactor(), umxHetCor(), umx_as_numeric(), umx_cont_2_quantiles(), umx_lower2full(), umx_make_MR_data(), umx_make_TwinData(), umx_make_fake_data(), umx_make_raw_from_cov(), umx_merge_randomized_columns(), umx_polychoric(), umx_polypairwise(), umx_polytriowise(), umx_read_lower(), umx_rename(), umx_reorder(), umx_score_scale(), umx_select_valid(), umx_stack(), umx_strings2numeric()

Examples

## Not run: 
fp = "~/Desktop/prolific_export_5f20c3e662e3b6407dcd37a5.csv"
df = prolific_read_demog(fp, sex = "Gender", age = "Age", df = df)
tmp = prolific_read_demog(fp, by.df = "PROLIFIC_PID", vars=c("Ethnicity.simplified"))

## End(Not run)

qm

Description

Quickmatrix function

Usage

qm(..., rowMarker = "|")

Arguments

Value

- matrix

See Also

Other Miscellaneous Utility Functions: install.OpenMx(), libs(), umx, umxLav2RAM(), umxModelNames(), umxRAM2Lav(), umxVersion(), umx_array_shift(), umx_find_object(), umx_lower.tri(), umx_msg(), umx_open_CRAN_page(), umx_pad(), umx_print()

Examples

# simple example
qm(0, 1 |
   2, NA)
## Not run: 
# clever example
M1 = M2 = diag(2)
qm(M1,c(4,5) | c(1,2),M2 | t(1:3))

## End(Not run)

Convert Radians to Degrees

Description

Just a helper to multiply radians by 180 and divide by \pi to get degrees.

note: R's trig functions, e.g. sin() use Radians for input!

There are 2\pi radians in a circle. 1 Rad = 180/\pi degrees = ~ 57.296 degrees.

Usage

rad2deg(rad)

Arguments

Value

• value in degrees

References

https://en.wikipedia.org/wiki/Radian

See Also

• deg2rad(), sin()

Other Miscellaneous Functions: deg2rad(), fin_JustifiedPE(), fin_NI(), fin_interest(), fin_percent(), fin_ticker(), fin_valuation(), umxBrownie()

Examples

rad2deg(pi) #180 degrees

Report coefficient alpha (reliability)

Description

Compute and report Coefficient alpha (extracted from Rcmdr to avoid its dependencies)

Usage

reliability(S)

Arguments

Value

None

References

- <https://cran.r-project.org/package=Rcmdr>

See Also

- [umx::print.reliability()],

Other Miscellaneous Stats Functions: FishersMethod(), SE_from_p(), geometric_mean(), harmonic_mean(), oddsratio(), umx, umxCov2cor(), umxHetCor(), umxParan(), umxWeightedAIC(), umx_apply(), umx_cor(), umx_means(), umx_r_test(), umx_round(), umx_scale(), umx_var()

Examples

# treat car data as items of a test
data(mtcars)
reliability(cov(mtcars))

Get residuals from an MxModel

Description

Return the residuals() from an OpenMx RAM model. You can format these (with digits), and suppress small values.

Usage

## S3 method for class 'MxModel'
residuals(object, digits = 2, suppress = NULL, reorder = NULL, ...)

Arguments

Value

• matrix of residuals

References

• https://tbates.github.io, https://github.com/tbates/umx

See Also

Other Reporting functions: RMSEA(), RMSEA.MxModel(), RMSEA.summary.mxmodel(), extractAIC.MxModel(), loadings(), loadings.MxModel(), tmx_show(), tmx_show.MxMatrix(), umxCI(), umxCI_boot(), umxConfint(), umxExpCov(), umxExpMeans(), umxFitIndices(), umxRotate()

Examples

## Not run: 
require(umx)
data(demoOneFactor)
manifests = names(demoOneFactor)

m1 = umxRAM("One Factor", data = demoOneFactor, type = "cov",
	umxPath("G", to = manifests),
	umxPath(var = manifests),
	umxPath(var = "G", fixedAt = 1.0)
)

# ===================================
# = Show the residuals of the model =
# ===================================
residuals(m1)
# |   |x1   |x2    |x3   |x4    |x5 |
# |:--|:----|:-----|:----|:-----|:--|
# |x1 |.    |.     |0.01 |.     |.  |
# |x2 |.    |.     |0.01 |-0.01 |.  |
# |x3 |0.01 |0.01  |.    |.     |.  |
# |x4 |.    |-0.01 |.    |.     |.  |
# |x5 |.    |.     |.    |.     |.  |
# [1] "nb: You can zoom in on bad values with, e.g. suppress = .01, which
#      will hide values smaller than this. Use digits = to round"

residuals(m1, digits = 3)
residuals(m1, digits = 3, suppress = .005)
# residuals are returned as an invisible object you can capture in a variable
a = residuals(m1); a

## End(Not run)

Graphical display of genotypic effects.

Description

tmx_genotypic_effect allows you to explore the concept of genotypic effect at a locus. With it, you can interactively explore the effects of allele frequency, additive variance, and dominance.

This function lets you explore the simplest two–allele system (B and b), with three possible genotypes, BB, Bb, and bb.

The point between the two homozygotes is m – the mean effect of the homozygous genotypes.

Parameter a is half the measured phenotypic difference between the homozygotes BB and bb. It corresponds to the additive effect of each additional B allele, relative to the bb phenotype.

Parameter d is the deviation of the heterozygote Bb phenotype from the homozygote mid-point m. It corresponds to the non-additive (dominance) effect of the B allele. The heterozygote phenotype may lie on either side of m and the sign of d will vary accordingly.

Old system from book ed 2:

Adapted from Mather and Jinks, 1977, p. 32). See book issue old-style nomenclature https://github.com/tbates/BGBook/issues/23

u = Frequency of the dominant allele (now = p). v = Frequency of the recessive allele (now = q).

m = midpoint between the two homozygotes d = half the difference between the two homozygote (now a)

h = deviation of the heterozygote from m (now = d)

New system:

u and v -> p and q

d and h -> a and d

See BGBook issue 23

Usage

tmx_genotypic_effect(p = 0.75, q = (1 - p), a = 0.5, d = 0, m = 0, show = TRUE)

Arguments

Value

• optional plot

References

• Neale, M. C. (2005). Quantitative Genetics. In Encyclopedia of Life Sciences. New York: John Wiley & Sons, Ltd. pdf

See Also

Other Teaching and testing Functions: tmx_is.identified(), umx

Examples

library(umx);


# =========================
# = Pure additivity: d= 0 =
# =========================
tmx_genotypic_effect(p = .5, a = 1, d = 0, m = 0, show = TRUE);

# =============================
# = Complete dominance: a=d=1 =
# =============================
tmx_genotypic_effect(p = .5, q =.5, a = 1, d = 1, m = 0, show = TRUE);

# ===========================
# = Over dominance: a< d =1 =
# ===========================
tmx_genotypic_effect(p = .5, q =.5, a =.5, d = 1, m = 0)

p = tmx_genotypic_effect(p = .5, q = .5, a = 1, d = .5, m = 0, show = TRUE); 
# p = p + ggplot2::geom_point() 
# p + ggplot2::geom_text(hjust = 0, nudge_x = 0.05, label= "x")
# ggsave(paste0(base, "c03_genotypic_effect_by_gene_dose.pdf"), width = 4.6, height = 4.6)

Test if a factor model is identified

Description

Test if a factor model is identified by establishing if the number of variables is equal too or grater than the number of model parameters. See also OpenMx::mxCheckIdentification() for checking actual models.

Usage

tmx_is.identified(nVariables, nFactors)

Arguments

Value

• Binary

References

• https://github.com/tbates/umx, https://tbates.github.io

See Also

• OpenMx::mxCheckIdentification()

Other Teaching and testing Functions: tmx_genotypic_effect(), umx

Examples

tmx_is.identified(nVariables = 2, nFactors = 1) # FALSE
tmx_is.identified(nVariables = 3, nFactors = 1) # TRUE
tmx_is.identified(nVariables = 4, nFactors = 2) # FALSE
tmx_is.identified(nVariables = 5, nFactors = 2) # TRUE

Show matrices of models in a easy-to-learn-from format.

Description

Show matrices of models in a easy-to-learn-from format.

Usage

tmx_show(
  x,
  what = c("values", "free", "labels", "nonzero_or_free"),
  show = c("free", "fixed", "all"),
  matrices = c("S", "A", "M"),
  digits = 2,
  report = c("html", "markdown"),
  na.print = "",
  zero.print = ".",
  html_font = NULL,
  style = c("paper", "material_dark", "classic", "classic_2", "minimal", "material"),
  bootstrap_options = c("hover", "bordered", "condensed", "responsive"),
  lightable_options = "striped"
)

Arguments

Value

None

See Also

Other Reporting functions: RMSEA(), RMSEA.MxModel(), RMSEA.summary.mxmodel(), extractAIC.MxModel(), loadings(), loadings.MxModel(), residuals.MxModel(), tmx_show.MxMatrix(), umxCI(), umxCI_boot(), umxConfint(), umxExpCov(), umxExpMeans(), umxFitIndices(), umxRotate()

Show matrices of models in a easy-to-learn-from format.

Description

Show matrices of models in a easy-to-learn-from format.

Usage

## S3 method for class 'MxMatrix'
tmx_show(
  x,
  what = c("values", "free", "labels", "nonzero_or_free"),
  show = c("free", "fixed", "all"),
  matrices = c("S", "A", "M"),
  digits = 2,
  report = c("html", "markdown"),
  na.print = "",
  zero.print = ".",
  html_font = NULL,
  style = c("paper", "material_dark", "classic", "classic_2", "minimal", "material"),
  bootstrap_options = c("hover", "bordered", "condensed", "responsive"),
  lightable_options = "striped"
)

Arguments

Value

None

See Also

Other Reporting functions: RMSEA(), RMSEA.MxModel(), RMSEA.summary.mxmodel(), extractAIC.MxModel(), loadings(), loadings.MxModel(), residuals.MxModel(), tmx_show(), umxCI(), umxCI_boot(), umxConfint(), umxExpCov(), umxExpMeans(), umxFitIndices(), umxRotate()

Show matrices of RAM models in a easy-to-learn-from format.

Description

A great way to learn about models is to look at the matrix contents. tmx_show is designed to do this in a way that makes it easy to process for users: The matrix contents are formatted as tables, and can even be displayed as tables in a web browser.

Usage

## S3 method for class 'MxModel'
tmx_show(
  x,
  what = c("values", "free", "labels", "nonzero_or_free"),
  show = c("free", "fixed", "all"),
  matrices = c("S", "A", "M"),
  digits = 2,
  report = c("html", "markdown"),
  na.print = "",
  zero.print = ".",
  html_font = NULL,
  style = c("paper", "material_dark", "classic", "classic_2", "minimal", "material"),
  bootstrap_options = c("hover", "bordered", "condensed", "responsive"),
  lightable_options = "striped"
)

Arguments

Details

The user can select which matrices to view, whether to show values, free, and/or labels, and the precision of rounding.

Value

None

References

• https://tbates.github.io

See Also

Other Teaching and Testing functions: umxDiagnose(), umxPower()

Examples

## Not run: 
require(umx)
data(demoOneFactor)
manifests = names(demoOneFactor)
m1 = umxRAM("tmx_sh", data = demoOneFactor, type = "cov",
	umxPath("G", to = manifests),
	umxPath(var = manifests),
	umxPath(var = "G", fixedAt = 1)
)

# =============================================
# = Show smart table on the web (the default) =
# =============================================
tmx_show(m1, report = "html")
tmx_show(m1, what = "free", matrices = "thresholds")
tmx_show(m1, zero.print = "-")

tmx_show(m1, report = "markdown")
tmx_show(m1, digits = 3, report = "markdown")
tmx_show(m1, matrices = "S", report = "markdown")
tmx_show(m1, what = "free"  , report = "markdown")
tmx_show(m1, what = "labels", report = "markdown")
tmx_show(m1, what = "free", matrices = "A", report= "markdown")

## End(Not run)

Deprecated. May already stop() code and ask to be updated. May be dropped entirely in future.

Description

xmuMakeThresholdsMatrices should be replaced with umxThresholdMatrix()

umxTryHard is deprecated: use umxRun() instead

stringToMxAlgebra is deprecated: please use umx_string_to_algebra() instead

genEpi_EvalQuote is deprecated: please use OpenMx::mxEvalByName() instead

umxReportCIs is deprecated: please use umxCI() instead

replace umxReportFit with umxSummary()

Replace umxGraph_RAM with plot()

Replace tryHard with OpenMx::mxTryHard()

Replace standardizeRAM with umx_standardize()

Arguments

References

• https://tbates.github.io, https://github.com/tbates/umx

Build and run a 2-group Cholesky ACE twin model (univariate or multivariate)

Description

Implementing a core task in twin modeling, umxACE models the genetic and environmental structure of one or more phenotypes (measured variables) using the Cholesky ACE model (Neale and Cardon, 1996).

Classical twin modeling uses the genetic and environmental differences among pairs of mono-zygotic (MZ) and di-zygotic (DZ) twins reared together.

umxACE implements a 2-group model to capture these data and represent the phenotypic variance as a sum of Additive genetic, unique environmental (E) and, optionally, either common or shared-environment (C) or non-additive genetic effects (D).

The following figure shows the ACE model for one variable "x" as a path diagram:

umxACE allows multivariate analyses, and this brings us to the Cholesky part of the model.

The Cholesky decomposition creates as many latent A (and C and E) latent variables as there are phenotypes, and, moving from left to right, decomposes the variance in each phenotype into successively restricted factors. The following figure shows the multivariate ACE model for three variables:

In this ACE model of three phenotypes, the expected variance-covariance matrix of the original data is the product of each lower Cholesky and its transform (i.e., A = a %*% t(a) summed for A+C+E.

This lower-triangle decomposition feature of the Cholesky yields a model which is certain to both account for all the variance (with some restrictions) in the data and be solvable.

This figure also contains the key to understanding how to modify models that umxACE produces using umxModify() to drop paths by label like "a_r1c1". nb: Read the "Matrices and Labels in ACE model" section in details below...

NOTE: Scroll down to details for how to use the function, a figure and multiple examples.

Usage

umxACE(
  name = "ACE",
  selDVs,
  selCovs = NULL,
  dzData = NULL,
  mzData = NULL,
  sep = NULL,
  data = NULL,
  zyg = "zygosity",
  type = c("Auto", "FIML", "cov", "cor", "WLS", "DWLS", "ULS"),
  numObsDZ = NULL,
  numObsMZ = NULL,
  boundDiag = 0,
  allContinuousMethod = c("cumulants", "marginals"),
  autoRun = getOption("umx_auto_run"),
  intervals = FALSE,
  tryHard = c("no", "yes", "ordinal", "search"),
  optimizer = NULL,
  residualizeContinuousVars = FALSE,
  nSib = 2,
  dzAr = 0.5,
  dzCr = 1,
  weightVar = NULL,
  equateMeans = TRUE,
  addStd = TRUE,
  addCI = TRUE
)

Arguments

Details

Covariates umxACE handles covariates by modelling them in the means. On the plus side, there is no distributional assumption for this method. A downside of this approach is that all covariates must be non-NA, thus dropping any rows where one or more covariates are missing. This can waste data. See also umx_residualize()).

Data Input The function flexibly accepts raw data, and also summary covariance data (in which case the user must also supple numbers of observations for the two input data sets).

The type parameter can select how you want the model data treated. "FIML" is the normal treatment. "cov" and "cor" will turn raw data into cor data for analysis, or check that you've provided cor data as input.

Types "WLS", "DWLS", and "ULS" will process raw data into WLS data of these types.

The default, "Auto" will treat data as the type they are provided as.

Ordinal Data In an important capability, the model transparently handles ordinal (binary or multi-level ordered factor data) inputs, and can handle mixtures of continuous, binary, and ordinal data in any combination. An experimental feature is under development to allow Tobit modeling.

The function also supports weighting of individual data rows. In this case, the model is estimated for each row individually, then each row likelihood is multiplied by its weight, and these weighted likelihoods summed to form the model-likelihood, which is to be minimized. This feature is used in the non-linear GxE model functions.

Additional features The umxACE function supports varying the DZ genetic association (defaulting to .5) to allow exploring assortative mating effects, as well as varying the DZ “C” factor from 1 (the default for modeling family-level effects shared 100% by twins in a pair), to .25 to model dominance effects.

Matrices and Labels in ACE model

Matrices 'a', 'c', and 'e' contain the path loadings of the Cholesky ACE factor model.

So, labels relevant to modifying the model are of the form ⁠"a_r1c1", "c_r1c1"⁠ etc.

Variables are in rows, and factors are in columns. So to drop the influence of factor 2 on variable 3, you would say:

m2 = umxModify(m1, update = "c_r3c2")

Less commonly-modified matrices are the mean matrix expMean. This has 1 row, and the columns are laid out for each variable for twin 1, followed by each variable for twin 2.

So, in a model where the means for twin 1 and twin 2 had been equated (set = to T1), you could make them independent again with this script:

m1$top$expMean$labels[1, 4:6] = c("expMean_r1c4", "expMean_r1c5", "expMean_r1c6")

note: Only one of C or D may be estimated simultaneously. This restriction reflects the lack of degrees of freedom to simultaneously model C and D with only MZ and DZ twin pairs (Eaves et al. 1978, p267).

Value

• OpenMx::mxModel() of subclass mxModel.ACE

References

• Eaves, L. J., Last, K. A., Young, P. A., & Martin, N. G. (1978). Model-fitting approaches to the analysis of human behaviour. Heredity, 41, 249-320. doi:10.1038/hdy.1978.101

See Also

• umxPlotACE(), umxSummaryACE(), power.ACE.test(), umxModify()

Other Twin Modeling Functions: power.ACE.test(), umx, umxACEcov(), umxACEv(), umxCP(), umxDiffMZ(), umxDiscTwin(), umxDoC(), umxDoCp(), umxGxE(), umxGxE_window(), umxGxEbiv(), umxIP(), umxMRDoC(), umxReduce(), umxReduceACE(), umxReduceGxE(), umxRotate.MxModelCP(), umxSexLim(), umxSimplex(), umxSummarizeTwinData(), umxSummaryACE(), umxSummaryACEv(), umxSummaryDoC(), umxSummaryGxEbiv(), umxSummarySexLim(), umxSummarySimplex(), umxTwinMaker()

Examples

require(umx)
# ============================
# = How heritable is height? =
# ============================

# 1. Height in meters has a tiny variance, and this makes optimising hard.
#    We'll scale it by 10x to make the Optimizer's task easier.
data(twinData) # ?twinData from Australian twins.
twinData[, c("ht1", "ht2")] = twinData[, c("ht1", "ht2")] * 10

# 2. Make mz & dz data.frames (no need to select variables: umx will do this)
mzData = twinData[twinData$zygosity %in% "MZFF", ]
dzData = twinData[twinData$zygosity %in% "DZFF", ]

# 3. Built & run the model, controlling for age in the means model
m1 = umxACE(selDVs = "ht", selCovs = "age", sep = "", dzData = dzData, mzData = mzData)

# sidebar: umxACE figures out variable names using sep: 
#    e.g. selVars = "wt" + sep= "_T" -> "wt_T1" "wt_T2"

umxSummary(m1, std = FALSE) # un-standardized

# tip 1: set report = "html" and umxSummary prints the table to your browser!
# tip 2: plot works for umx: Get a figure of the model and parameters
# plot(m1) # Also, look at the options for ?plot.MxModel.

# ===========================================
# = Test ADE, AE, CE, E, and generate table =
# ===========================================

umxReduce(m1, report="html", silent= TRUE)

# ============================
# = Model, with 2 covariates =
# ============================

# Create another covariate: cohort
twinData$cohort1 = twinData$cohort2 =twinData$part
mzData = twinData[twinData$zygosity %in% "MZFF", ]
dzData = twinData[twinData$zygosity %in% "DZFF", ]

# 1. def var approach
m2 = umxACE(selDVs = "ht", selCovs = c("age", "cohort"), sep = "", dzData = dzData, mzData = mzData)

# 2. Residualized approach: remove height variance accounted-for by age.
FFdata = twinData[twinData$zygosity %in% c("MZFF", "DZFF"), ]
FFdata = umx_residualize("ht", "age", suffixes = 1:2, data = FFdata)
mzData = FFdata[FFdata$zygosity %in% "MZFF", ]
dzData = FFdata[FFdata$zygosity %in% "DZFF", ]
m3 = umxACE(selDVs = "ht", sep = "", dzData = dzData, mzData = mzData)

# =============================================================
# = ADE: Evidence for dominance ? (DZ correlation set to .25) =
# =============================================================
m2 = umxACE(selDVs = "ht", sep = "", dzData = dzData, mzData = mzData, dzCr = .25)
umxCompare(m2, m1) # ADE is better
umxSummary(m2, comparison = m1) 
# nb: Although summary is smart enough to print d, the underlying 
#     matrices are still called a, c & e.

# tip: try umxReduce(m1) to automatically build and compare ACE, ADE, AE, CE
# including conditional probabilities!

# ===================================================
# = WLS example using diagonal weight least squares =
# ===================================================

m3 = umxACE(selDVs = "ht", sep = "", dzData = dzData, mzData = mzData, 
	type = "DWLS", allContinuousMethod='marginals'
)


# ==============================
# = Univariate model of weight =
# ==============================

# Things to note:

# 1. Weight has a large variance, and this makes solution finding very hard.
# Here, we residualize the data for age, which also scales weight and height.

data(twinData)
tmp = umx_residualize(c("wt", "ht"), cov = "age", suffixes= c(1, 2), data = twinData)
mzData = tmp[tmp$zygosity %in% "MZFF", ]
dzData = tmp[tmp$zygosity %in% "DZFF", ]

# tip: You might also want transform variables
# tmp = twinData$wt1[!is.na(twinData$wt1)]
# car::powerTransform(tmp, family="bcPower"); hist(tmp^-0.6848438)
# twinData$wt1 = twinData$wt1^-0.6848438
# twinData$wt2 = twinData$wt2^-0.6848438

# 4. note: the default boundDiag = 0 lower-bounds a, c, and e at 0.
#    Prevents mirror-solutions. If not desired: set boundDiag = NULL.

m2 = umxACE(selDVs = "wt", dzData = dzData, mzData = mzData, sep = "", boundDiag = NULL)

# A short cut (which is even shorter for "_T" twin data with "MZ"/"DZ" data in zygosity column is:
m1 = umxACE(selDVs = "wt", sep = "", data = twinData,
	dzData = c("DZMM", "DZFF", "DZOS"), mzData = c("MZMM", "MZFF"))
# |   |   a1|c1 |   e1|
# |:--|----:|:--|----:|
# |wt | 0.93|.  | 0.38|

# tip: umx_make_twin_data_nice() will make data into this nice format for you!

# ======================
# = MODEL MODIFICATION =
# ======================
# We can modify this model, e.g. test shared environment. 
# Set comparison to modify, and show effect in one step.

m2 = umxModify(m1, update = "c_r1c1", name = "no_C", comparison = TRUE)
#*tip* call umxModify(m1) with no parameters, and it will print the labels available to fix!
# nb: You can see parameters of any model with parameters(m1)

# =========================================================
# = Well done! Now you can make modify twin models in umx =
# =========================================================

# =====================================
# = Bivariate height and weight model =
# =====================================
data(twinData)
# We'll scale height (ht1 and ht2) and weight
twinData = umx_scale_wide_twin_data(data = twinData, varsToScale = c("ht", "wt"), sep = "")
mzData = twinData[twinData$zygosity %in% c("MZFF", "MZMM"),]
dzData = twinData[twinData$zygosity %in% c("DZFF", "DZMM", "DZOS"), ]
m1 = umxACE(selDVs = c("ht", "wt"), sep = '', dzData = dzData, mzData = mzData)
umxSummary(m1)

# ===================
# = Ordinal example =
# ===================

# Prep data
require(umx)
data(twinData)
# Cut BMI column to form ordinal obesity variables
obLevels = c('normal', 'overweight', 'obese')
cuts = quantile(twinData[, "bmi1"], probs = c(.5, .2), 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 umxFactors
ordDVs = c("obese1", "obese2")
twinData[, ordDVs] = umxFactor(twinData[, ordDVs])

mzData = twinData[twinData$zygosity %in% "MZFF", ]
dzData = twinData[twinData$zygosity %in% "DZFF", ]

# Model and summary!
m1 = umxACE(selDVs = "obese", dzData = dzData, mzData = mzData, sep = '')

# And controlling age (otherwise manifests appearance as latent C)
m1 = umxACE(selDVs = "obese", selCov= "age", dzData = dzData, mzData = mzData, sep = '')
# umxSummary(m1)

# ============================================
# = Bivariate continuous and ordinal example =
# ============================================
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, .2), 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])
mzData = twinData[twinData$zygosity %in% "MZFF",] 
dzData = twinData[twinData$zygosity %in% "DZFF",]
mzData = mzData[1:80,] # just top 80 so example runs in a couple of secs
dzData = dzData[1:80,]
m1 = umxACE(selDVs= c("wt","obese"), dzData= dzData, mzData= mzData, sep='')

# And controlling age
m1 = umxACE(selDVs = c("wt","obese"), selCov= "age", dzData = dzData, mzData = mzData, sep = '')

# =======================================
# = Mixed continuous and binary example =
# =======================================
require(umx)
data(twinData)
twinData= umx_scale_wide_twin_data(data= twinData,varsToScale= "wt", sep="")
# Cut to form category of 20% obese subjects
# and make into mxFactors (ensure ordered is TRUE, and require levels)
obLevels   = c('normal', 'obese')
cuts       = quantile(twinData[, "bmi1"], probs = .2, 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) 
ordDVs = c("obese1", "obese2")
twinData[, ordDVs] = umxFactor(twinData[, ordDVs])

selDVs = c("wt", "obese")
mzData = twinData[twinData$zygosity %in% "MZFF",]
dzData = twinData[twinData$zygosity %in% "DZFF",]
m1 = umxACE(selDVs = selDVs, dzData = dzData, mzData = mzData, sep = '')
umxSummary(m1)

# ==============
# = Two binary =
# ==============
require(umx)
data(twinData)
htLevels   = c('short', 'tall')
obLevels   = c('normal', 'obese')
cuts       = quantile(twinData[, "bmi1"], probs = .2, 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) 
ordDVs = c("obese1", "obese2")
twinData[, ordDVs] = umxFactor(twinData[, ordDVs])

twinData$short1 = cut(twinData$ht1, breaks=c(-Inf,1.6,Inf), labels=htLevels) 
twinData$short2 = cut(twinData$ht2, breaks=c(-Inf,1.6,Inf), labels=htLevels) 
ordDVs = c("short1", "short2")
twinData[, ordDVs] = umxFactor(twinData[, ordDVs])

mzData = twinData[twinData$zygosity %in% "MZFF",]
dzData = twinData[twinData$zygosity %in% "DZFF",]
m1 = umxACE(selDVs = c("short", "obese"), dzData = dzData, mzData = mzData, sep = '')

# ===================================
# Example with covariance data only =
# ===================================

require(umx)
data(twinData)
twinData= umx_scale_wide_twin_data(data=twinData, varsToScale= "wt", sep="")
selDVs = c("wt1", "wt2")
mz = cov(twinData[twinData$zygosity %in%  "MZFF", selDVs], use = "complete")
dz = cov(twinData[twinData$zygosity %in%  "DZFF", selDVs], use = "complete")
m1 = umxACE(selDVs=selDVs, dzData=dz, mzData=mz, numObsDZ=569, numObsMZ=351)
umxSummary(m1)
plot(m1)

Run a Cholesky with covariates that are random (in the expected covariance matrix)

Description

Often, researchers include covariates in 2-group Cholesky umxACE() twin models. The umxACEcov 'random' option models the covariates in the expected covariance matrix, thus allowing all data to be preserved. The downside is that this method has a strong assumption of multivariate normality. Covariates like age, which are perfectly correlated in twins cannot be used. Covariates like sex, which are ordinal, violate the normality assumption. Binary and ordinal covariates like sex also violate the normality assumption. Which is most of the use cases :-(.

Usage

umxACEcov(
  name = "ACEcov",
  selDVs,
  selCovs,
  dzData,
  mzData,
  sep = NULL,
  type = c("Auto", "FIML", "cov", "cor", "WLS", "DWLS", "ULS"),
  allContinuousMethod = c("cumulants", "marginals"),
  dzAr = 0.5,
  dzCr = 1,
  addStd = TRUE,
  addCI = TRUE,
  boundDiag = 0,
  equateMeans = TRUE,
  bVector = FALSE,
  autoRun = getOption("umx_auto_run"),
  tryHard = c("no", "yes", "ordinal", "search"),
  optimizer = NULL
)

Arguments

Details

The following figure shows how the ACE model with random covariates appears as a path diagram:

Value

• OpenMx::mxModel() of subclass mxModel.ACEcov

References

• Neale, M. C., & Martin, N. G. (1989). The effects of age, sex, and genotype on self-report drunkenness following a challenge dose of alcohol. Behavior Genetics, 19, 63-78. doi:10.1007/BF01065884.

• Schwabe, I., Boomsma, D. I., Zeeuw, E. L., & Berg, S. M. (2015). A New Approach to Handle Missing Covariate Data in Twin Research : With an Application to Educational Achievement Data. Behavior Genetics, 46, 583-95. doi:10.1007/s10519-015-9771-1.

See Also

Other Twin Modeling Functions: power.ACE.test(), umx, umxACE(), umxACEv(), umxCP(), umxDiffMZ(), umxDiscTwin(), umxDoC(), umxDoCp(), umxGxE(), umxGxE_window(), umxGxEbiv(), umxIP(), umxMRDoC(), umxReduce(), umxReduceACE(), umxReduceGxE(), umxRotate.MxModelCP(), umxSexLim(), umxSimplex(), umxSummarizeTwinData(), umxSummaryACE(), umxSummaryACEv(), umxSummaryDoC(), umxSummaryGxEbiv(), umxSummarySexLim(), umxSummarySimplex(), umxTwinMaker()

Examples

## Not run: 
# ============================================
# = BMI, can't use Age as a random covariate =
# ============================================
require(umx)
data(twinData)
# Replicate age to age1 & age2
twinData$age1 = twinData$age2 = twinData$age
mzData = subset(twinData, zygosity == "MZFF")
dzData = subset(twinData, zygosity == "DZFF")

# =====================================================================
# = Trying to use identical var (like age) as a random cov is ILLEGAL =
# =====================================================================
m1 = umxACEcov(selDVs = "bmi", selCovs = "age", dzData = dzData, mzData = mzData, sep = "")

# ========================================================
# = Use an lm-based age-residualisation approach instead =
# ========================================================

resid_data = umx_residualize("bmi", "age", suffixes = 1:2, twinData)
mzData = subset(resid_data, zygosity == "MZFF")
dzData = subset(resid_data, zygosity == "DZFF")
m2     = umxACE("resid", selDVs = "bmi", dzData = dzData, mzData = mzData, sep = "")

# Univariate BMI without covariate of age for comparison
mzData = subset(twinData, zygosity == "MZFF")
dzData = subset(twinData, zygosity == "DZFF")
m3 = umxACE("raw_bmi", selDVs = "bmi", dzData = dzData, mzData = mzData, sep = "")

# ===========================================================================
# = A bivariate example (need a dataset with a VIABLE COVARIATE to do this) =
# ===========================================================================
selDVs  = "wt" # Set the DVs
selCovs = "ht" # Set the COV
selVars = umx_paste_names(selDVs, covNames = selCovs, sep = "", sep = 1:2)
mzData = subset(twinData, zygosity == "MZFF")
dzData = subset(twinData, zygosity == "DZFF")
m1 = umxACEcov(selDVs = selDVs, selCovs = selCovs,
   dzData = dzData, mzData = mzData, sep = "", autoRun = TRUE
)

## End(Not run)

Build and run 2-group uni- or multi-variate ACE models based on VARIANCE (not paths).

Description

A common task in twin modeling involves using the genetic and environmental differences between large numbers of pairs of mono-zygotic (MZ) and di-zygotic (DZ) twins reared together to model the genetic and environmental structure of one, or, typically, several phenotypes. umxACEv directly estimates variance components (rather than paths, which are then squared to produce variance and therefore cannot be negative). It offers better power, correct Type I error and un-biased estimates (with no zero-bound for the variances) as a saturated model. (Verhulst et al, 2019).

The ACE variance-based model decomposes phenotypic variance into additive genetic (A), unique environmental (E) and, optionally, either common environment (shared-environment, C) or non-additive genetic effects (D). Scroll down to details for how to use the function, a figure and multiple examples.

The following figure shows the A components of a trivariate ACEv model:

NOTE: This function does not use the Cholesky decomposition. Instead it directly models variance. This ensures unbiased type-I error rates. It means that occasionally estimates of variance may be negative. This should be used as an occasion to inspect you model choices and data. umxACEv can be used as a base model to validate the ACE Cholesky model, a core model in behavior genetics (Neale and Cardon, 1992).

Usage

umxACEv(
  name = "ACEv",
  selDVs,
  selCovs = NULL,
  sep = NULL,
  dzData,
  mzData,
  dzAr = 0.5,
  dzCr = 1,
  type = c("Auto", "FIML", "cov", "cor", "WLS", "DWLS", "ULS"),
  allContinuousMethod = c("cumulants", "marginals"),
  data = NULL,
  zyg = "zygosity",
  weightVar = NULL,
  numObsDZ = NULL,
  numObsMZ = NULL,
  addStd = TRUE,
  addCI = TRUE,
  boundDiag = NULL,
  equateMeans = TRUE,
  bVector = FALSE,
  autoRun = getOption("umx_auto_run"),
  tryHard = c("no", "yes", "ordinal", "search"),
  optimizer = NULL,
  nSib = 2
)

Arguments

Details

Data Input The function flexibly accepts raw data, and also summary covariance data (in which case the user must also supple numbers of observations for the two input data sets).

Ordinal Data In an important capability, the model transparently handles ordinal (binary or multi-level ordered factor data) inputs, and can handle mixtures of continuous, binary, and ordinal data in any combination.

The function also supports weighting of individual data rows. In this case, the model is estimated for each row individually, then each row likelihood is multiplied by its weight, and these weighted likelihoods summed to form the model-likelihood, which is to be minimized. This feature is used in the non-linear GxE model functions.

Additional features The umxACEv function supports varying the DZ genetic association (defaulting to .5) to allow exploring assortative mating effects, as well as varying the DZ “C” factor from 1 (the default for modeling family-level effects shared 100% by twins in a pair), to .25 to model dominance effects.

note: Only one of C or D may be estimated simultaneously. This restriction reflects the lack of degrees of freedom to simultaneously model C and D with only MZ and DZ twin pairs (Eaves et al. 1978 p267).

Value

• OpenMx::mxModel() subclass mxModelACEv

References

• Verhulst, B., Prom-Wormley, E., Keller, M., Medland, S., & Neale, M. C. (2019). Type I Error Rates and Parameter Bias in Multivariate Behavioral Genetic Models. Behav Genet, 49, 99-111. doi:10.1007/s10519-018-9942-y Eaves, L. J., Last, K. A., Young, P. A., & Martin, N. G. (1978). Model-fitting approaches to the analysis of human behaviour. Heredity, 41, 249-320. doi:10.1038/hdy.1978.101

See Also

Other Twin Modeling Functions: power.ACE.test(), umx, umxACE(), umxACEcov(), umxCP(), umxDiffMZ(), umxDiscTwin(), umxDoC(), umxDoCp(), umxGxE(), umxGxE_window(), umxGxEbiv(), umxIP(), umxMRDoC(), umxReduce(), umxReduceACE(), umxReduceGxE(), umxRotate.MxModelCP(), umxSexLim(), umxSimplex(), umxSummarizeTwinData(), umxSummaryACE(), umxSummaryACEv(), umxSummaryDoC(), umxSummaryGxEbiv(), umxSummarySexLim(), umxSummarySimplex(), umxTwinMaker()

Examples

## Not run: 

# ==============================
# = Univariate model of weight =
# ==============================
require(umx)
data(twinData) # ?twinData from Australian twins.

# Things to note: ACE model of weight will return a NEGATIVE variance in C.
#  This is exactly why we have ACEv! It suggests we need a different model
#  In this case: ADE.
# Other things to note:
# 1. umxACEv can figure out variable names: provide "sep", and selVars. 
#    Function generates: "wt" -> "wt1" "wt2"
# 2. umxACEv picks the variables it needs from the data.

mzData = twinData[twinData$zygosity %in% "MZFF", ]
dzData = twinData[twinData$zygosity %in% "DZFF", ]
m1 = umxACEv(selDVs = "wt", sep = "", dzData = dzData, mzData = mzData)

# A short cut (which is even shorter for "_T" twin data with "MZ"/"DZ" data in zygosity column is:
m1 = umxACEv(selDVs = "wt", sep = "", dzData = "MZFF", mzData = "DZFF", data = twinData)
# ========================================================
# = Evidence for dominance ? (DZ correlation set to .25) =
# ========================================================
m2 = umxACEv("ADE", selDVs = "wt", sep = "", dzData = dzData, mzData = mzData, dzCr = .25)
# note: the underlying matrices are still called A, C, and E.
# I catch this in the summary table, so columns are labeled A, D, and E.
# However, currently, the plot will say A, C, E.

# We can modify this model, dropping dominance component (still called C), 
# and see a comparison:
m3 = umxModify(m2, update = "C_r1c1", comparison = TRUE, name="AE")
# =========================================================
# = Well done! Now you can make modify twin models in umx =
# =========================================================

# ============================
# = How heritable is height? =
# ============================
# 
# Note: Height has a small variance. umx can typically picks good starts,
#    but scaling is advisable.
# 
require(umx)
# Load data and rescale height to cm (var in m too small)
data(twinData) # ?twinData from Australian twins.
twinData[,c("ht1", "ht2")]= twinData[,c("ht1", "ht2")]*100

mzData = twinData[twinData$zygosity %in% "MZFF", ]
dzData = twinData[twinData$zygosity %in% "DZFF", ]
m1 = umxACEv(selDVs = "ht", sep = "", dzData = dzData, mzData = mzData)

umxSummary(m1, std = FALSE) # unstandardized
plot(m1)

# tip: with report = "html", umxSummary can print the table to your browser!
# tip: You can turn off auto-plot with umx_set_auto_plot(FALSE)

# ========================================================
# = Evidence for dominance ? (DZ correlation set to .25) =
# ========================================================
m2 = umxACEv("ADE", selDVs = "ht", dzCr = .25, sep="", dzData = dzData, mzData = mzData)
umxCompare(m2, m1) # Is ADE better?
umxSummary(m2, comparison = m1) # nb: though this is ADE, matrices are still called A,C,E

# We can modify this model, dropping shared environment, and see a comparison:
m3 = umxModify(m2, update = "C_r1c1", comparison = TRUE, name = "AE")

# =====================================
# = Bivariate height and weight model =
# =====================================

data(twinData)
twinData[,c("ht1", "ht2")]= twinData[,c("ht1", "ht2")]*100
mzData = twinData[twinData$zygosity %in% c("MZFF", "MZMM"), ]
dzData = twinData[twinData$zygosity %in% c("DZFF", "DZMM", "DZOS"), ]
m1 = umxACEv(selDVs = c("ht", "wt"), sep = '', dzData = dzData, mzData = mzData)

# ===================
# = Ordinal example =
# ===================
require(umx)
data(twinData)

# Cut bmi column to form ordinal obesity variables
cutPoints = quantile(twinData[, "bmi1"], probs = c(.5, .2), na.rm = TRUE)
obesityLevels = c('normal', 'overweight', 'obese')
twinData$obese1 = cut(twinData$bmi1, breaks = c(-Inf, cutPoints, Inf), labels = obesityLevels) 
twinData$obese2 = cut(twinData$bmi2, breaks = c(-Inf, cutPoints, Inf), labels = obesityLevels) 

# Make the ordinal variables into mxFactors (ensure ordered is TRUE, and require levels)
twinData[, c("obese1", "obese2")] = umxFactor(twinData[, c("obese1", "obese2")])
mzData = twinData[twinData$zygosity %in% "MZFF", ]
dzData = twinData[twinData$zygosity %in% "DZFF", ]
m2 = umxACEv(selDVs = "obese", dzData = dzData, mzData = mzData, sep = '')

# FYI: Show mz, dz, and t1 and t2 have the same levels!
str(mzData)

# ============================================
# = Bivariate continuous and ordinal example =
# ============================================
data(twinData)
# Cut bmi column to form ordinal obesity variables
ordDVs = c("obese1", "obese2")
obesityLevels = c('normal', 'overweight', 'obese')
cutPoints = quantile(twinData[, "bmi1"], probs = c(.5, .2), na.rm = TRUE)
twinData$obese1 = cut(twinData$bmi1, breaks = c(-Inf, cutPoints, Inf), labels = obesityLevels) 
twinData$obese2 = cut(twinData$bmi2, breaks = c(-Inf, cutPoints, Inf), labels = obesityLevels) 

# Make the ordinal variables into ordered mxFactors
twinData[, ordDVs] = umxFactor(twinData[, ordDVs])

# umxACEv can trim out unused variables on its own
mzData = twinData[twinData$zygosity %in% "MZFF", ]
dzData = twinData[twinData$zygosity %in% "DZFF", ]

m1 = umxACEv(selDVs = c("wt", "obese"), dzData = dzData, mzData = mzData, sep = '')
plot(m1)

# =======================================
# = Mixed continuous and binary example =
# =======================================
require(umx)
data(twinData)
# Cut to form category of 20% obese subjects
# and make into mxFactors (ensure ordered is TRUE, and require levels)
cutPoints = quantile(twinData[, "bmi1"], probs = .2, na.rm = TRUE)
obesityLevels = c('normal', 'obese')
twinData$obese1 = cut(twinData$bmi1, breaks = c(-Inf, cutPoints, Inf), labels = obesityLevels) 
twinData$obese2 = cut(twinData$bmi2, breaks = c(-Inf, cutPoints, Inf), labels = obesityLevels) 
ordDVs = c("obese1", "obese2")
twinData[, ordDVs] = umxFactor(twinData[, ordDVs])

selDVs = c("wt", "obese")
mzData = twinData[twinData$zygosity %in% "MZFF", ]
dzData = twinData[twinData$zygosity %in% "DZFF", ]
m1 = umxACEv(selDVs = selDVs, dzData = dzData, mzData = mzData, sep = '')
umxSummary(m1)

# ===================================
# Example with covariance data only =
# ===================================

require(umx)
data(twinData)
selDVs = c("wt")
mz = cov(twinData[twinData$zygosity %in% "MZFF", tvars(selDVs, "")], use = "complete")
dz = cov(twinData[twinData$zygosity %in% "DZFF", tvars(selDVs, "")], use = "complete")
m1 = umxACEv(selDVs = selDVs, sep= "", dzData = dz, mzData= mz, numObsDZ= 569, numObsMZ= 351)
umxSummary(m1, std = FALSE)

## End(Not run)

Creates nicely formatted journal-style summaries of models, p-values, data-frames and much more.

Description

umxAPA creates APA-style reports from a range of statistical models, or to summarize data. I wrote it to suit me.

Nice alternatives include jtools::summ.

Example functionality includes:

1. Given an stats::lm() model, umxAPA will return a formatted effect, including 95% CI. e.g.: umxAPA(lm(mpg~wt, data=mtcars), "wt") yields: \beta = -5.34 [-6.48, -4.20], p < 0.001. here "wt" restricts the output to just the named effect.

2. umxAPA also supports t.test(), stats::glm(), cor.test(), and others as I need them.

3. Get a CI from obj=beta and se=se : umxAPA(-0.30, .03) returns \beta = -0.3 [-0.36, -0.24]

4. Back out an SE from \beta and CI: umxAPA(-0.030, c(-0.073, 0.013)) returns \beta = -0.03, se = 0.02

5. Given only a number as obj, will be treated as a p-value, and returned in APA format.

6. Given a dataframe, umxAPA will return a table of correlations with means and SDs in the last row. e.g.: ⁠umxAPA(mtcars[,c("cyl", "wt", "mpg", )]⁠ yields:

   	cyl	wt	mpg
   cyl	1	0.78	-0.85
   wt	0.78	1	-0.87
   mpg	-0.85	-0.87	1
   mean_sd	6.19 (1.79)	3.22 (0.98)	20.09 (6.03)

Usage

umxAPA(
  obj = .Last.value,
  se = NULL,
  p = NULL,
  std = FALSE,
  digits = 2,
  use = "complete",
  min = 0.001,
  addComparison = NA,
  report = c("markdown", "html", "none", "expression"),
  lower = TRUE,
  test = c("Chisq", "LRT", "Rao", "F", "Cp"),
  SEs = TRUE,
  means = TRUE,
  suffix = "",
  cols = NA
)

Arguments

Value

• string

References

• https://stats.oarc.ucla.edu/r/dae/logit-regression/

See Also

SE_from_p()

Other Reporting Functions: umx, umxFactorScores(), umxGetLatents(), umxGetManifests(), umxGetModel(), umxGetParameters(), umxParameters(), umx_aggregate(), umx_time()

Examples

# ========================================
# = Report lm (regression/anova) results =
# ========================================
umxAPA(lm(mpg ~ wt + disp, mtcars)) # Report all parameters
umxAPA(lm(mpg ~ wt + disp, mtcars), "wt") # Just effect of weight
umxAPA(lm(mpg ~ wt + disp, mtcars), std = TRUE) # Standardize model!

###############
# GLM example #
###############

df = mtcars
df$mpg_thresh = 0
df$mpg_thresh[df$mpg > 16] = 1
m1 = glm(mpg_thresh ~ wt + gear,data = df, family = binomial)
umxAPA(m1)

###############
# A t-Test    #
###############

umxAPA(t.test(x = 1:10, y = c(7:20)))
umxAPA(t.test(extra ~ group, data = sleep))

# ======================================================
# = Summarize DATA FRAME: Correlations + Means and SDs =
# ======================================================
umxAPA(mtcars[,1:3])
umxAPA(mtcars[,1:3], digits = 3)
umxAPA(mtcars[,1:3], lower = FALSE)
## Not run: 
umxAPA(mtcars[,1:3], report = "html")

## End(Not run)

# ==========================================
# = CONFIDENCE INTERVAL from effect and se =
# ==========================================
umxAPA(.4, .3) # parameter 2 interpreted as SE

# Input beta and CI, and back out the SE
umxAPA(-0.030, c(-0.073, 0.013), digits = 3)

# ====================
# = Format a p-value =
# ====================
umxAPA(.0182613)   #   0.02
umxAPA(.00018261) # < 0.001
umxAPA(.00018261, addComparison = FALSE) # 0.001

# ========================
# = Report a correlation =
# ========================
data(twinData)
tmp = subset(twinData, zygosity %in% c("MZFF", "MZMM"))
m1 = cor.test(~ wt1 + wt2, data = tmp)
umxAPA(m1)

A simple wrapper for mxAlgebra with name as the first parameter for more readable compact code.

Description

umxAlgebra is a wrapper for mxAlgebra which has the name parameter first in order.

Usage

umxAlgebra(
  name = NA,
  expression,
  dimnames = NA,
  ...,
  joinKey = as.character(NA),
  joinModel = as.character(NA),
  verbose = 0L,
  initial = matrix(as.numeric(NA), 1, 1),
  recompute = c("always", "onDemand"),
  fixed = "deprecated_use_recompute"
)

Arguments

Value

• OpenMx::mxAlgebra()

See Also

• umxMatrix()

Other Advanced Model Building Functions: umx, umxFixAll(), umxJiggle(), umxRun(), umxThresholdMatrix(), umxUnexplainedCausalNexus(), xmuLabel(), xmuValues()

Examples

## Not run: 
A = umxMatrix("A", "Full", nrow = 3, ncol = 3, values=2)
B = umxAlgebra("B", A)
C = umxAlgebra(A + B, name = "C")
D = umxAlgebra(sin(C), name = "D")
m1 = mxRun(mxModel("AlgebraExample", A, B, C, D ))
mxEval(D, m1)
	
x = umxAlgebra("circ", expression = 2 * pi)
class(x$formula)
x = mxAlgebra(name = "circ", 2 * pi)
class(x$formula) # "call"

## End(Not run)

A recipe Easter-egg for umx

Description

How to cook steak.

Usage

umxBrownie()

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.

See Also

• OpenMx::omxBrownie()

Other Miscellaneous Functions: deg2rad(), fin_JustifiedPE(), fin_NI(), fin_interest(), fin_percent(), fin_ticker(), fin_valuation(), rad2deg()

Examples

umxBrownie()

Add (and, optionally, run) confidence intervals to a structural model.

Description

umxCI adds OpenMx::mxCI() calls for requested (default all) parameters in a model, runs these CIs if necessary, and reports them in a neat summary.

Usage

umxCI(
  model = NULL,
  which = c("ALL", NA, "list of your making"),
  remove = FALSE,
  run = c("no", "yes", "if necessary", "show"),
  interval = 0.95,
  type = c("both", "lower", "upper"),
  regex = NULL,
  showErrorCodes = TRUE
)

Arguments

Details

umxCI also reports if any problems were encountered. The codes are standard OpenMx errors and warnings

• 1: The final iterate satisfies the optimality conditions to the accuracy requested, but the sequence of iterates has not yet converged. NPSOL was terminated because no further improvement could be made in the merit function (Mx status GREEN)

• 2: The linear constraints and bounds could not be satisfied. The problem has no feasible solution.

• 3: The nonlinear constraints and bounds could not be satisfied. The problem may have no feasible solution.

• 4: The major iteration limit was reached (Mx status BLUE).

• 6: The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)

• 7: The function derivatives returned by funcon or funobj appear to be incorrect.

• 9: An input parameter was invalid.

If run = "no", the function simply adds the CI requests, but returns the model without running them.

Value

• OpenMx::mxModel()

References

• https://github.com/tbates/umx

See Also

• stats::confint(), umxConfint(), umxCI(), umxModify()

Other Reporting functions: RMSEA(), RMSEA.MxModel(), RMSEA.summary.mxmodel(), extractAIC.MxModel(), loadings(), loadings.MxModel(), residuals.MxModel(), tmx_show(), tmx_show.MxMatrix(), umxCI_boot(), umxConfint(), umxExpCov(), umxExpMeans(), umxFitIndices(), umxRotate()

Examples

## Not run: 
require(umx)
data(demoOneFactor)
manifests = names(demoOneFactor)

m1 = umxRAM("One Factor", data = demoOneFactor, type = "cov",
	umxPath("G", to = manifests),
	umxPath(var = manifests),
	umxPath(var = "G", fixedAt = 1)
)
m1$intervals # none yet - empty list()
m1 = umxCI(m1)
m1$intervals # $G_to_x1...
m1 = umxCI(m1, remove = TRUE) # remove CIs from the model and return it
m1$intervals # none again

# Add CIs by name
parameters(m1, patt="_with_")
m1 = umxCI(m1, which = "x1_with_x1")
m1 = umxCI(m1, which = c("x1_with_x1", "x2_with_x2"))
m1 = umxCI(m1, regex = "x1_with_", run= "yes")
#           lbound estimate ubound lbound Code ubound Code
# x1_with_x1  0.036    0.041  0.047           0           0

# ========================
# = A twin model example =
# ========================
data(twinData) 
mzData = subset(twinData, zygosity == "MZFF")
dzData = subset(twinData, zygosity == "DZFF")
m1 = umxACE(selDVs = c("bmi1","bmi2"), dzData = dzData, mzData = mzData)
umxCI(m1, run = "show") # show what will be requested
umxCI(m1, run = "yes") # actually compute the CIs
# Don't force update of CIs, but if they were just added, then calculate them
umxCI(m1, run = "if necessary")
m1 = umxCI(m1, remove = TRUE) # remove them all
m1$intervals # none!
# Show what parameters are available to get CIs on
umxParameters(m1) 
# Request a CI by label:
m1 = umxCI(m1, which = "a_r1c1", run = "yes")

## End(Not run)

umxCI_boot

Description

Compute boot-strapped Confidence Intervals for parameters in an OpenMx::mxModel() The function creates a sampling distribution for parameters by repeatedly drawing samples with replacement from your data and then computing the statistic for each redrawn sample.

Usage

umxCI_boot(
  model,
  rawData = NULL,
  type = c("par.expected", "par.observed", "empirical"),
  std = TRUE,
  rep = 1000,
  conf = 95,
  dat = FALSE,
  digits = 3
)

Arguments

Value

• expected covariance matrix

See Also

• umxExpMeans(), umxExpCov()

Other Reporting functions: RMSEA(), RMSEA.MxModel(), RMSEA.summary.mxmodel(), extractAIC.MxModel(), loadings(), loadings.MxModel(), residuals.MxModel(), tmx_show(), tmx_show.MxMatrix(), umxCI(), umxConfint(), umxExpCov(), umxExpMeans(), umxFitIndices(), umxRotate()

Examples

## Not run: 
require(umx)
data(demoOneFactor)
manifests = names(demoOneFactor)

m1 = umxRAM("One Factor", data = demoOneFactor, type = "cov",
	umxPath("G", to = manifests),
	umxPath(var = manifests),
	umxPath(var = "G", fixedAt = 1.0)
)

umxCI_boot(m1, type = "par.expected")

## End(Not run)

Runs cross-lagged panel models

Description

One way of assessing causal relationships is by introducing time into the analyses. umxCLPM implements three cross-lagged panel models (CLPM) from the literature. The first is the classic CLPM from Heise (1969), the second is the CLPM from Hamaker et al. (2015), and the third is the CLPM from STARTS (1995). You simply pass the number of waves and the data set along with the model you wish to run.

Sketch mode is available; if you pass column names to data, a model object is returned for manipulation later.

Usage

umxCLPM(
  waves,
  name = NULL,
  model = c("Hamaker2015", "Heise1969", "STARTS1995", "IV_RI_CLPM"),
  data = NULL,
  counts = NULL,
  summary = !umx_set_silent(silent = TRUE),
  autoRun = getOption("umx_auto_run"),
  tryHard = c("no", "yes", "ordinal", "search"),
  verbose = FALSE,
  batteries = c("scale", "ordinaloptim"),
  std = FALSE,
  ivs = NULL,
  defn = NULL,
  defto = NULL,
  type = c("Auto", "FIML", "cov", "cor", "WLS", "DWLS", "ULS"),
  allContinuousMethod = c("cumulants", "marginals")
)

Arguments

Value

An OpenMx model object.

References

Kenny, D.A., & Zautra, A. (1995). The trait-state-error model for multiwave data. Journal of Consulting and Clinical Psychology, 63, 52–59. doi:10.1037/0022-006X.63.1.52 Hamaker E.L., Kuiper R.M., & Grasman R. (2015). A critique of the cross-lagged panel model. Psychological Methods, 20, 102–116. doi:10.1037/a0038889 Heise D. R. (1970). Causal inference from panel data. Sociological Methodology, 2, 3–27. doi:10.2307/270780

Examples

## Not run: 

# ================
# = 1. Load Data =
# ================
data(docData)
dt <- docData[2:9]

# ============================
# = 2. Make a CLPM model     =
# ============================
hamaker <- umxCLPM(waves = 4, name = "mymodel", model = "Hamaker2015", data = dt)

## End(Not run)

umxCP: Build and run a Common Pathway twin model

Description

Make a 2-group Common Pathway twin model.

The common-pathway model (aka "psychometric model" (McArdle and Goldsmith, 1990) provides a powerful tool for theory-based testing of genetic and environmental differences. It proposes that A, C, and E components act on a latent substrate (organ, mental mechanism etc.) and this is manifested in the measured phenotypes.

umxCP supports this with pairs of mono-zygotic (MZ) and di-zygotic (DZ) twins reared together to model the genetic and environmental structure of multiple phenotypes (measured behaviors).

Common-pathway path diagram:

As can be seen, each phenotype also by default has A, C, and E influences specific to that phenotype.

Features include the ability to include more than one common pathway, to use ordinal data.

note: The function umx_set_optimization_options() allows users to see and set mvnRelEps and mvnMaxPointsA mvnRelEps defaults to .005. For ordinal models, you might find that '0.01' works better.

Usage

umxCP(
  name = "CP",
  selDVs,
  selCovs = NULL,
  dzData = NULL,
  mzData = NULL,
  sep = NULL,
  nFac = 1,
  type = c("Auto", "FIML", "cov", "cor", "WLS", "DWLS", "ULS"),
  data = NULL,
  zyg = "zygosity",
  allContinuousMethod = c("cumulants", "marginals"),
  correlatedACE = FALSE,
  dzAr = 0.5,
  dzCr = 1,
  autoRun = getOption("umx_auto_run"),
  tryHard = c("yes", "no", "ordinal", "search"),
  optimizer = NULL,
  equateMeans = TRUE,
  weightVar = NULL,
  bVector = FALSE,
  boundDiag = 0,
  addStd = TRUE,
  addCI = TRUE,
  numObsDZ = NULL,
  numObsMZ = NULL,
  freeLowerA = FALSE,
  freeLowerC = FALSE,
  freeLowerE = FALSE,
  correlatedA = "deprecated"
)

Arguments

Details

Like the umxACE() model, the CP model decomposes phenotypic variance into additive genetic (A), unique environmental (E) and, optionally, either common or shared-environment (C) or non-additive genetic effects (D).

Unlike the Cholesky, these factors do not act directly on the phenotype. Instead latent A, C, and E influences impact on one or more latent factors which in turn account for variance in the phenotypes (see Figure).

Data Input Currently, the umxCP function accepts only raw data. This may change in future versions.

Ordinal Data

In an important capability, the model transparently handles ordinal (binary or multi-level ordered factor data) inputs, and can handle mixtures of continuous, binary, and ordinal data in any combination.

Additional features

The umxCP function supports varying the DZ genetic association (defaulting to .5) to allow exploring assortative mating effects, as well as varying the DZ “C” factor from 1 (the default for modeling family-level effects shared 100% by twins in a pair), to .25 to model dominance effects.

Matrices and Labels in CP model

A good way to see which matrices are used in umxCP is to run an example model and plot it.

All the shared matrices are in the model "top".

Matrices top$as, top$cs, and top$es contain the path loadings specific to each variable on their diagonals.

So, to see the 'as' values, labels, or free states, you can say:

m1$top$as$values

m1$top$as$free

m1$top$as$labels

Labels relevant to modifying the specific loadings take the form "as_r1c1", "as_r2c2" etc.

The common-pathway loadings on the factors are in matrices top$a_cp, top$c_cp, top$e_cp.

The common factors themselves are in the matrix top$cp_loadings (an nVar * 1 matrix)

Less commonly-modified matrices are the mean matrix expMean. This has 1 row, and the columns are laid out for each variable for twin 1, followed by each variable for twin 2. So, in a model where the means for twin 1 and twin 2 had been equated (set = to T1), you could make them independent again with this line:

m1$top$expMean$labels[1,4:6] = c("expMean_r1c4", "expMean_r1c5", "expMean_r1c6")

For a deep-dive, see xmu_make_TwinSuperModel()

Value

• OpenMx::mxModel()

References

• Martin, N. G., & Eaves, L. J. (1977). The Genetical Analysis of Covariance Structure. Heredity, 38, 79-95.

• Kendler, K. S., Heath, A. C., Martin, N. G., & Eaves, L. J. (1987). Symptoms of anxiety and symptoms of depression. Same genes, different environments? Archives of General Psychiatry, 44, 451-457. doi:10.1001/archpsyc.1987.01800170073010.

• McArdle, J. J., & Goldsmith, H. H. (1990). Alternative common factor models for multivariate biometric analyses. Behavior Genetics, 20, 569-608. doi:10.1007/BF01065873.

• https://github.com/tbates/umx

See Also

• umxSummaryCP(), umxPlotCP(). See umxRotate.MxModelCP() to rotate the factor loadings of a umxCP() model. See umxACE() for more examples of twin modeling. plot() and umxSummary() work for all twin models, e.g., umxIP(), umxCP(), umxGxE(), and umxACE().

Other Twin Modeling Functions: power.ACE.test(), umx, umxACE(), umxACEcov(), umxACEv(), umxDiffMZ(), umxDiscTwin(), umxDoC(), umxDoCp(), umxGxE(), umxGxE_window(), umxGxEbiv(), umxIP(), umxMRDoC(), umxReduce(), umxReduceACE(), umxReduceGxE(), umxRotate.MxModelCP(), umxSexLim(), umxSimplex(), umxSummarizeTwinData(), umxSummaryACE(), umxSummaryACEv(), umxSummaryDoC(), umxSummaryGxEbiv(), umxSummarySexLim(), umxSummarySimplex(), umxTwinMaker()

Examples

## Not run: 
# ========================================================
# = Run a 3-factor Common pathway twin model of 6 traits =
# ========================================================
require(umx)
data(GFF)
mzData = subset(GFF, zyg_2grp == "MZ")
dzData = subset(GFF, zyg_2grp == "DZ")
selDVs = c("gff", "fc", "qol", "hap", "sat", "AD") 
m1 = umxCP(selDVs = selDVs, sep = "_T", nFac = 3, tryHard = "yes",
		dzData = dzData, mzData = mzData)

# Shortcut using "data ="
selDVs = c("gff", "fc", "qol", "hap", "sat", "AD") 
m1 = umxCP(selDVs= selDVs, nFac= 3, data=GFF, zyg="zyg_2grp")

# ===================
# = Do it using WLS =
# ===================
m2 = umxCP("new", selDVs = selDVs, sep = "_T", nFac = 3, optimizer = "SLSQP",
		dzData = dzData, mzData = mzData, tryHard = "ordinal", 
	type= "DWLS", allContinuousMethod='marginals'
)

# =================================================
# = Find and test dropping of shared environment  =
# =================================================
# Show all labels for C parameters 
umxParameters(m1, patt = "^c")
# Test dropping the 9 specific and common-factor C paths
m2 = umxModify(m1, regex = "(cs_.*$)|(c_cp_)", name = "dropC", comp = TRUE)
umxSummaryCP(m2, comparison = m1, file = NA)
umxCompare(m1, m2)

# =======================================
# = Mixed continuous and binary example =
# =======================================
data(GFF)
# Cut to form umxFactor 20% depressed  DEP
cutPoints = quantile(GFF[, "AD_T1"], probs = .2, na.rm = TRUE)
ADLevels  = c('normal', 'depressed')
GFF$DEP_T1 = cut(GFF$AD_T1, breaks = c(-Inf, cutPoints, Inf), labels = ADLevels) 
GFF$DEP_T2 = cut(GFF$AD_T2, breaks = c(-Inf, cutPoints, Inf), labels = ADLevels) 
ordDVs = c("DEP_T1", "DEP_T2")
GFF[, ordDVs] = umxFactor(GFF[, ordDVs])

selDVs = c("gff","fc","qol","hap","sat","DEP") 
mzData = subset(GFF, zyg_2grp == "MZ")
dzData = subset(GFF, zyg_2grp == "DZ")

# umx_set_optimizer("NPSOL")
# umx_set_optimization_options("mvnRelEps", .01)
m1 = umxCP(selDVs = selDVs, sep = "_T", nFac = 3, dzData = dzData, mzData = mzData)
m2 = umxModify(m1, regex = "(cs_r[3-5]|c_cp_r[12])", name = "dropC", comp= TRUE)

# Do it using WLS
m3 = umxCP(selDVs = selDVs, sep = "_T", nFac = 3, dzData = dzData, mzData = mzData,
		tryHard = "ordinal", type= "DWLS")
# TODO umxCPL fix WLS here
# label at row 1 and column 1 of matrix 'top.binLabels'' in model 'CP3fac' : object 'Vtot'

# ==============================
# = Correlated factors example =
# ==============================
# ====================
# = DON'T USE THIS!!! =
# ====================
data(GFF)
mzData = subset(GFF, zyg_2grp == "MZ")
dzData = subset(GFF, zyg_2grp == "DZ")
selDVs = c("gff", "fc", "qol", "hap", "sat", "AD")
m1 = umxCP("base_model", selDVs = selDVs, sep = "_T", correlatedACE = TRUE, 
	 dzData = dzData, mzData = mzData, nFac = 3, tryHard = "yes")

# What are the ace covariance labels? (two ways to get)
umx_lower.tri(m1$top$a_cp$labels)
parameters(m1, patt = "[ace]_cp")

# 1. Now allow a1 and a2 to correlate
m2=umxModify(m1,regex="a_cp_r2c1",name="a2_a1_cov",free=TRUE,tryHard="yes")
umxCompare(m2, m1)

# 2. Drop all (a|c|e) correlations from a model
tmp= namez(umx_lower.tri(m2$top$a_cp$labels), "a_cp", replace= "[ace]_cp")
m3 = umxModify(m2, regex= tmp, comparison = TRUE)

## End(Not run) # end dontrun

Print a comparison table of one or more OpenMx::mxModel()s, formatted nicely.

Description

umxCompare compares two or more OpenMx::mxModel()s. It has several nice features:

1. It supports direct control of rounding, and reports p-values rounded to APA style.

2. It reports the table in your preferred format (default is markdown, options include latex).

3. Table columns are arranged to make for easy comparison for readers.

4. report = 'inline', will provide an English sentence suitable for a paper.

5. report = "html" opens a web table in your browser to paste into a word processor.

Note: If you leave comparison blank, it will just give fit info for the base model

Usage

umxCompare(
  base = NULL,
  comparison = NULL,
  all = TRUE,
  digits = 3,
  report = c("markdown", "html", "inline"),
  compareWeightedAIC = FALSE,
  silent = FALSE,
  file = "tmp.html"
)

Arguments

References

• https://github.com/tbates/umx

See Also

• umxSummary(), umxRAM(),umxCompare()

Other Model Summary and Comparison: umx, umxEquate(), umxMI(), umxReduce(), umxSetParameters(), umxSummary()

Examples

## Not run: 
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)
)

m2 = umxModify(m1, update = "G_to_x2", name = "drop_path_2_x2")
umxCompare(m1, m2)
umxCompare(m1, m2, report = "inline") # Add English-sentence descriptions
umxCompare(m1, m2, report = "html") # Open table in browser

# Two comparison models
m3 = umxModify(m2, update = "G_to_x3", name = "drop_path_2_x2_and_3")

umxCompare(m1, c(m2, m3))
umxCompare(m1, c(m2, m3), compareWeightedAIC = TRUE)
umxCompare(c(m1, m2), c(m2, m3), all = TRUE)

manifests = names(demoOneFactor)
m1 = umxRAM("WLS", data = demoOneFactor, type = "DWLS",
	umxPath("G", to = manifests),
	umxPath(var = manifests),
	umxPath(var = "G", fixedAt = 1)
)

m2 = umxModify(m1, update = "G_to_x2", name = "drop_path_2_x2")
umxCompare(m1, m2)
umxCompare(m1, m2, report = "inline") # Add English-sentence descriptions
umxCompare(m1, m2, report = "html") # Open table in browser

## End(Not run)

Get confidence intervals from a umx model

Description

Implements confidence interval function for umx models.

Usage

umxConfint(
  object,
  parm = c("existing", "all", "or one or more labels", "smart"),
  wipeExistingRequests = TRUE,
  level = 0.95,
  run = FALSE,
  showErrorCodes = FALSE,
  optimizer = c("SLSQP", "NPSOL", "CSOLNP", "current")
)

Arguments

Details

Note: By default, requesting new CIs wipes the existing ones. To keep these, set wipeExistingRequests = FALSE.

Because CIs can take time to run, by default only already-computed CIs will be reported. To run new CIs, set run = TRUE .

Note: OpenMx defines a confint function which will return SE-based CIs.

If parm is empty, and run = FALSE, a message will alert you to set run = TRUE.

Value

• OpenMx::mxModel()

References

• https://github.com/tbates/umx

See Also

• stats::confint(), OpenMx::mxSE(), umxCI(), OpenMx::mxCI()

Other Reporting functions: RMSEA(), RMSEA.MxModel(), RMSEA.summary.mxmodel(), extractAIC.MxModel(), loadings(), loadings.MxModel(), residuals.MxModel(), tmx_show(), tmx_show.MxMatrix(), umxCI(), umxCI_boot(), umxExpCov(), umxExpMeans(), umxFitIndices(), umxRotate()

Examples

## Not run: 
require(umx)
data(demoOneFactor)

manifests = names(demoOneFactor)
m1 = umxRAM("OneFactor", data = demoOneFactor, type = "cov",
	umxPath(from = "G", to = manifests),
	umxPath(var = manifests),
	umxPath(var = "G", fixedAt = 1)
)

m1 = umxConfint(m1, run = TRUE) # There are no existing CI requests...

# Add a CI request for "G_to_x1", run, and report. Save with this CI computed
m2 = umxConfint(m1, parm = "G_to_x1", run = TRUE) 

# Just print out any existing CIs
umxConfint(m2)

# CI requests added for free matrix parameters. User prompted to set run = TRUE
m3 = umxConfint(m1, "all")

# Run the requested CIs
m3 = umxConfint(m3, run = TRUE) 

# Run CIs for free one-headed (asymmetric) paths in RAM model. 
#   note: Deletes other existing requests,
tmp = umxConfint(m1, parm = "A", run = TRUE)

# Wipe existing CIs, add G_to_x1
tmp = umxConfint(m1, parm = "G_to_x1", run = TRUE, wipeExistingRequests = TRUE) 

# For some twin models, a "smart" mode is implemented
# note: only implemented for umxCP so far
m2 =  umxConfint(m1, "smart")

## End(Not run)

Convert a covariance matrix into a correlation matrix

Description

A version of cov2cor() that forces upper and lower triangles to be identical (rather than nearly identical)

Usage

umxCov2cor(x)

Arguments

Value

• A correlation matrix

References

• https://github.com/tbates/umx

See Also

cov2cor()

Other Miscellaneous Stats Functions: FishersMethod(), SE_from_p(), geometric_mean(), harmonic_mean(), oddsratio(), reliability(), umx, umxHetCor(), umxParan(), umxWeightedAIC(), umx_apply(), umx_cor(), umx_means(), umx_r_test(), umx_round(), umx_scale(), umx_var()

Examples

umxCov2cor(cov(mtcars[,1:5]))

Diagnose problems in a model - not working!

Description

The goal of this function WILL BE (not currently functional) to diagnose problems in a model and return suggestions to the user. It is a work in progress, and of no use as yet.

Usage

umxDiagnose(model, tryHard = FALSE, diagonalizeExpCov = FALSE)

Arguments

Details

Best diagnostics are:

1. Observed data variances and means

2. Expected variances and means

3. Difference of these?

Try * diagonalizeExpCov diagonal * umx_is_ordered()

Tricky, but reporting variances and standardized thresholds is ideal. Guidance is to start with unit variances and thresholds within +/- 2 SD of the mean. Like %p option in Classic Mx.

Value

• helpful messages and perhaps a modified model

References

• https://tbates.github.io, https://github.com/tbates/umx

See Also

Other Teaching and Testing functions: tmx_show.MxModel(), umxPower()

Examples

## Not run: 
require(umx)
data(demoOneFactor)
manifests = names(demoOneFactor)

m1 = umxRAM("OneFactor", data = demoOneFactor, type= "cov",
	umxPath("G", to = manifests),
	umxPath(var = manifests),
	umxPath(var = "G", fixedAt = 1)
)
m1 = mxRun(m1)
umxSummary(m1, std = TRUE)
umxDiagnose(m1)

## End(Not run)

MZ differences method for testing evidence for causality.

Description

umxDiffMZ implements the simple twin1-twin2 based correlation method, e.g. De Moor (2008), in which MZ differences on a variable x asserted to be causal of an outcome variable y are tested for association with differences on y. The logic of the design is shown below:

Usage

umxDiffMZ(
  x,
  y,
  data,
  sep = "_T",
  mzZygs = c("MZFF", "MZMM"),
  zyg = "zygosity",
  labxy = c(-1.2, 1.8),
  xylim = c(NA, NA),
  digits = 2
)

Arguments

Details

Example output is shown below, with the fitted line and fit inscribed. The plot is just a ggplot graph that is returned and can be edited and formatted.

For a more sophisticated linear mixed model approach, see umxDiscTwin().

Value

• Graph for decorating

References

• De Moor, M. H., Boomsma, D. I., Stubbe, J. H., Willemsen, G., & de Geus, E. J. (2008). Testing causality in the association between regular exercise and symptoms of anxiety and depression. Archives of General Psychiatry, 65(8), 897-905. doi:10.1001/archpsyc.65.8.897.

See Also

• umxDoC(), umxDiscTwin(), umxMR()

Other Twin Modeling Functions: power.ACE.test(), umx, umxACE(), umxACEcov(), umxACEv(), umxCP(), umxDiscTwin(), umxDoC(), umxDoCp(), umxGxE(), umxGxE_window(), umxGxEbiv(), umxIP(), umxMRDoC(), umxReduce(), umxReduceACE(), umxReduceGxE(), umxRotate.MxModelCP(), umxSexLim(), umxSimplex(), umxSummarizeTwinData(), umxSummaryACE(), umxSummaryACEv(), umxSummaryDoC(), umxSummaryGxEbiv(), umxSummarySexLim(), umxSummarySimplex(), umxTwinMaker()

Examples

data(twinData)
umxDiffMZ(x="ht", y="wt", labxy = c(-.5, 3), data = twinData, sep = "")
umxDiffMZ(x="ht", y="wt", xylim = c( -2, 2), data = twinData, sep = "")

Intra-pair association in MZ, DZ twin models. (ALPHA quality!)

Description

Testing causal claims is often difficult due to an inability to experimentally randomize traits and situations. A combination of control data and data from twins discordant for the putative causal trait can falsify causal hypotheses.

umxDiscTwin uses nlme::nlme() to compute the beta for x in y ~ x in models either a) Only controlling non-independence, and b) MZ and DZ subsample models in which the family level of the predictor y is also controlled.

If x is causal, then the effect size of x on y is expected to be equally large in all three samples. If the population association reflects confounded genes or shared environments, then the association in MZ twins will reduce to zero/non-significance.

The function uses the nlme::lme() function to compute the effect of the presumed causal variable on the outcome, controlling, for mid-family score and with random means model using familyID. e.g.:

mzModel = lme(fixed = y ~ x + FamMeanX, random = ~ 1+FamMeanX|FAMID, data = umx_scale(MZ), na.action = "na.omit")

Example output from umxDiscTwin

Usage

umxDiscTwin(
  x,
  y,
  data,
  mzZygs = c("MZFF", "MZMM"),
  dzZygs = c("DZFF", "DZMM", "DZOS"),
  FAMID = "FAMID",
  out = c("table", "plot", "model"),
  use = "complete.obs",
  sep = "_T"
)

Arguments

Value

• table of results

References

• Begg, M. D., & Parides, M. K. (2003). Separation of individual-level and cluster-level covariate effects in regression analysis of correlated data. Stat Med, 22(16), 2591-2602. doi:10.1002/sim.1524

• Bergen, S. E., Gardner, C. O., Aggen, S. H., & Kendler, K. S. (2008). Socioeconomic status and social support following illicit drug use: causal pathways or common liability? Twin Res Hum Genet, 11, 266-274. doi:10.1375/twin.11.3.266

• McGue, M., Osler, M., & Christensen, K. (2010). Causal Inference and Observational Research: The Utility of Twins. Perspectives on Psychological Science, 5, 546-556. doi:10.1177/1745691610383511

See Also

• umxDoC(), umxDiffMZ(), umxMR()

Other Twin Modeling Functions: power.ACE.test(), umx, umxACE(), umxACEcov(), umxACEv(), umxCP(), umxDiffMZ(), umxDoC(), umxDoCp(), umxGxE(), umxGxE_window(), umxGxEbiv(), umxIP(), umxMRDoC(), umxReduce(), umxReduceACE(), umxReduceGxE(), umxRotate.MxModelCP(), umxSexLim(), umxSimplex(), umxSummarizeTwinData(), umxSummaryACE(), umxSummaryACEv(), umxSummaryDoC(), umxSummaryGxEbiv(), umxSummarySexLim(), umxSummarySimplex(), umxTwinMaker()

Examples

## Not run: 
data(twinData)
# add to test must set FAMID umxDiscTwin(x = "ht", y = "wt", data = twinData, sep="")
tmp = umxDiscTwin(x = "ht", y = "wt", data = twinData, sep="", FAMID = "fam")
print(tmp, digits = 3)

## End(Not run)

Build and run a 2-group Direction of Causation twin models.

Description

Testing causal claims is often difficult due to an inability to conduct experimental randomization of traits and situations to people. When twins are available, even when measured on a single occasion, the pattern of cross-twin cross-trait correlations can (given distinguishable modes of inheritance for the two traits) falsify causal hypotheses.

umxDoC implements a 2-group model to form latent variables for each of two traits, and allows testing whether trait 1 causes trait 2, vice-versa, or even reciprocal causation.

Using latent variables instead of a manifest measure for testing causation, avoids the bias created by differences in measurement error in which the more reliable measure appears to "cause" the less reliable one (Gillespie and Martin, 2005).

The following figure shows how the DoC model appears as a path diagram (for two latent variables X and Y, each with three indicators). Note: For pedagogical reasons, only the model for 1 twin is shown, and only one DoC pathway drawn.

Usage

umxDoC(
  name = "DoC",
  var1Indicators,
  var2Indicators,
  mzData = NULL,
  dzData = NULL,
  sep = "_T",
  causal = TRUE,
  autoRun = getOption("umx_auto_run"),
  intervals = FALSE,
  tryHard = c("no", "yes", "ordinal", "search"),
  optimizer = NULL,
  data = NULL,
  zyg = "zygosity"
)

Arguments

Value

• OpenMx::mxModel() of subclass MxModelDoC

References

• N.A. Gillespie and N.G. Martin (2005). Direction of Causation Models. In Encyclopedia of Statistics in Behavioral Science, 1. 496–499. Eds. Brian S. Everitt & David C. Howell.

• McGue, M., Osler, M., & Christensen, K. (2010). Causal Inference and Observational Research: The Utility of Twins. Perspectives on Psychological Science, 5, 546-556. doi:10.1177/1745691610383511

• Rasmussen, S. H. R., Ludeke, S., & Hjelmborg, J. V. B. (2019). A major limitation of the direction of causation model: non-shared environmental confounding. Twin Res Hum Genet, 22, 1-13. doi:10.1017/thg.2018.67

See Also

• umxDiscTwin()

Other Twin Modeling Functions: power.ACE.test(), umx, umxACE(), umxACEcov(), umxACEv(), umxCP(), umxDiffMZ(), umxDiscTwin(), umxDoCp(), umxGxE(), umxGxE_window(), umxGxEbiv(), umxIP(), umxMRDoC(), umxReduce(), umxReduceACE(), umxReduceGxE(), umxRotate.MxModelCP(), umxSexLim(), umxSimplex(), umxSummarizeTwinData(), umxSummaryACE(), umxSummaryACEv(), umxSummaryDoC(), umxSummaryGxEbiv(), umxSummarySexLim(), umxSummarySimplex(), umxTwinMaker()

Examples

## Not run: 

# ========================
# = Does Rain cause Mud? =
# ========================

# ================
# = 1. Load Data =
# ================
data(docData)
docData = umx_scale_wide_twin_data(c(var1, var2), docData, sep= "_T")
mzData  = subset(docData, zygosity %in% c("MZFF", "MZMM"))
dzData  = subset(docData, zygosity %in% c("DZFF", "DZMM"))

# =======================================
# = 2. Define manifests for var 1 and 2 =
# =======================================
var1 = paste0("varA", 1:3)
var2 = paste0("varB", 1:3)

# =======================================================
# = 3. Make the non-causal (Cholesky) and causal models =
# =======================================================
Chol = umxDoC(var1= var1, var2= var2, mzData= mzData, dzData= dzData, causal= FALSE)
# nb: DoC initially has causal paths fixed @0
DoC  = umxDoC(var1= var1, var2= var2, mzData= mzData, dzData= dzData, causal= TRUE)
a2b   = umxModify(DoC, "a2b", free = TRUE, name = "a2b"); summary(a2b)
b2a   = umxModify(DoC, "b2a", free = TRUE, name = "b2a"); summary(b2a)
Recip = umxModify(DoC, c("a2b", "b2a"), free = TRUE, name = "Recip"); summary(Recip)

# Compare fits
umxCompare(Chol, c(a2b, b2a, Recip))

# ==========================================
# = Alternative call with data in one file =
# ==========================================
data(docData)
docData = umx_scale_wide_twin_data(c(var1, var2), docData, sep= "_T")
DoC = umxDoC(var1= paste0("varA", 1:3), var2= paste0("varB", 1:3),
	  mzData= c("MZFF", "MZMM"), dzData= c("DZFF", "DZMM"), data = docData
)

## End(Not run)

Make a direction of causation model based on umxPath statements

Description

Makes a direction of causation model with umxPath() statements

Usage

umxDoCp(
  var1Indicators,
  var2Indicators,
  mzData = NULL,
  dzData = NULL,
  sep = "_T",
  causal = TRUE,
  name = "DoC",
  autoRun = getOption("umx_auto_run"),
  intervals = FALSE,
  tryHard = c("no", "yes", "ordinal", "search"),
  optimizer = NULL
)

Arguments

Details

See also umxDoC()

Value

• [A direction of causation model with umxPath() statements.

See Also

• umxDoC()

Other Twin Modeling Functions: power.ACE.test(), umx, umxACE(), umxACEcov(), umxACEv(), umxCP(), umxDiffMZ(), umxDiscTwin(), umxDoC(), umxGxE(), umxGxE_window(), umxGxEbiv(), umxIP(), umxMRDoC(), umxReduce(), umxReduceACE(), umxReduceGxE(), umxRotate.MxModelCP(), umxSexLim(), umxSimplex(), umxSummarizeTwinData(), umxSummaryACE(), umxSummaryACEv(), umxSummaryDoC(), umxSummaryGxEbiv(), umxSummarySexLim(), umxSummarySimplex(), umxTwinMaker()

Examples

## Not run: 
# ================
# = 1. Load Data =
# ================
data(docData)
var1 = paste0("varA", 1:3)
var2 = paste0("varB", 1:3)
tmp = umx_scale_wide_twin_data(varsToScale= c(var1, var2), sep= "_T", data= docData)
mzData = subset(docData, zygosity %in% c("MZFF", "MZMM"))
dzData = subset(docData, zygosity %in% c("DZFF", "DZMM"))
m1 = umxDoCp(var1, var2, mzData= mzData, dzData= dzData, sep = "_T", causal= TRUE)


## End(Not run)

FIML-based Exploratory Factor Analysis (EFA)

Description

Perform full-information maximum-likelihood factor analysis on a data matrix.

Usage

umxEFA(
  x = NULL,
  factors = NULL,
  data = NULL,
  scores = c("none", "ML", "WeightedML", "Regression"),
  minManifests = NA,
  rotation = c("varimax", "promax", "none"),
  return = c("model", "loadings"),
  report = c("markdown", "html"),
  summary = FALSE,
  name = "efa",
  digits = 2,
  tryHard = c("no", "yes", "ordinal", "search"),
  n.obs = NULL,
  covmat = NULL
)

Arguments

Details

As in factanal(), you need only specify the number of factors and offer up some manifest data, e.g:

umxEFA(factors = 2, data = mtcars)

Equivalently, you can also give a list of factor names:

umxEFA(factors = c("g", "v"), data = mtcars)

The factor model is implemented as a structural equation model, e.g.

You can request scores from the model. Unlike factanal, these can cope with missing data.

You can also rotate the factors using any rotation function.

In an EFA, all items may load on all factors.

Should work with rotations provided in libs("GPArotation") and libs("psych"), e.g.,

Orthogonal: "varimax", "quartimax", "bentlerT", "equamax", "varimin", "geominT" and "bifactor" Oblique: "Promax", "promax", "oblimin", "simplimax", "bentlerQ", "geominQ", "biquartimin" and "cluster"

For identification we need m² degrees of freedom. We get m(m+1)/2 from fixing factor variances to 1 and covariances to 0. We get another m(m-1)/2 degrees of freedom by fixing the upper-right hand corner of the factor loadings component of the A matrix at 0.

To aid optimization, manifest residual variances are lbounded at 0.

EFA reports standardized loadings: to do this, we scale the data.

note: Bear in mind that factor scores are indeterminate (can be rotated to an infinity of equivalent solutions).

Thanks to @ConorDolan for code implementing the rotation matrix and other suggestions!

Value

• EFA OpenMx::mxModel()

References

• https://github.com/tbates/umx,

Hendrickson, A. E. and White, P. O. (1964). Promax: a quick method for rotation to orthogonal oblique structure. British Journal of Statistical Psychology, 17, 65–70. doi:10.1111/j.2044-8317.1964.tb00244.x.

Kaiser, H. F. (1958). The varimax criterion for analytic rotation in factor analysis. Psychometrika, 23, 187–200. doi:10.1007/BF02289233.

See Also

• factanal(), OpenMx::mxFactorScores()

Other Super-easy helpers: umx, umxTwoStage()

Examples

## Not run: 
myVars = c("mpg", "disp", "hp", "wt", "qsec")
m1 = umxEFA(mtcars[, myVars], factors =   2, rotation = "promax")
# By default, returns the model
umx_is_MxModel(m1) # TRUE
# The loadings are stashed in the model:
loadings(m1)

# Formula interface in umxEFA
m2 = umxEFA(~ mpg + disp + hp + wt + qsec, factors = 2, rotation = "promax", data = mtcars)
loadings(m2)

# base-R factanal Formula interface for comparison
m2 = factanal(~ mpg + disp + hp + wt + qsec, factors = 2, rotation = "promax", data = mtcars)
loadings(m2)

# Return the loadings object
x = umxEFA(mtcars[, myVars], factors = 2, return = "loadings")
names(x) # "loadings" "rotmat"

# scores requested, so these will be returned
x = umxEFA(name = "score", factors = "g", data = mtcars[, myVars], scores= "Regression")
head(x)
#       g
# 1  -0.48059346
# 2  -0.42354000
# 3  -0.87078110

m1 = umxEFA(myVars, factors = 2, data = mtcars, rotation = "promax")
m1 = umxEFA(name = "named", factors = "g", data = mtcars[, myVars])
m1 = umxEFA(name = "by_number", factors = 2, rotation = "promax", data = mtcars[, myVars])


## End(Not run)

umxEquate: Equate two or more paths

Description

In addition to dropping or adding parameters, a second common task in modeling is to equate parameters. umx provides a convenience function to equate parameters by setting one or more parameters (the "slave" set) equal to one or more "master" parameters. These parameters are picked out via their labels, and setting two or more parameters to have the same value is accomplished by setting the slave(s) to have the same label(s) as the master parameters, thus constraining them to take the same value during model fitting.

Usage

umxEquate(
  model,
  a,
  b,
  newlabels = NULL,
  free = c(TRUE, FALSE, NA),
  verbose = FALSE,
  name = NULL,
  autoRun = FALSE,
  tryHard = c("no", "yes", "ordinal", "search"),
  comparison = TRUE,
  master = NULL,
  slave = NULL
)

Arguments

Details

note: In addition to using this method to equating parameters, you can also equate one parameter to another by setting its label to the "square bracket" address of the master, e.g. "a[r,c]".

Tip: To find labels of free parameters use umxGetParameters() with free = TRUE

Tip: To find labels by name, use the regex parameter of umxGetParameters()

Value

• OpenMx::mxModel()

References

• https://github.com/tbates/umx

See Also

umxModify(), umxCompare()

Other Model Summary and Comparison: umx, umxCompare(), umxMI(), umxReduce(), umxSetParameters(), umxSummary()

Examples

## Not run: 
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)
)
# By default, umxEquate just equates master and slave labels: doesn't run model
m2 = umxEquate(m1, a = "G_to_x1", b = "G_to_x2", name = "Eq x1 x2 loadings")

# Set autoRun = TRUE and comparison = TRUE to run and output a comparison
m2 = umxEquate(m1, autoRun = TRUE, comparison = TRUE, name = "Eq_x1_x2",
	     a = "G_to_x1", b = "G_to_x2"
)

# rename the equated paths
m2 = umxEquate(m1, autoRun = TRUE, comparison = TRUE, name = "Eq_x1_x2",
	     a = "G_to_x1", b = "G_to_x2", newlabels = c("equated")
)
parameters(m2)

## End(Not run)

Example code from Twin Research and Human Genetics Paper on umx

Description

This is the example code used in our Twin Research and Human Genetics Paper on umx

Usage

umxExamples()

References

• Bates, T. C., Neale, M. C., & Maes, H. H. (2019). umx: A library for Structural Equation and Twin Modelling in R. Twin Research and Human Genetics, 22, 27-41. doi:10.1017/thg.2019.2.

See Also

• umx()

Examples

## Not run: 

# ==========================================================================
# = Example code from Twin Research and Human Genetics Paper on umx(model) =
# ==========================================================================

# Installing umx can be done using the R-code:
install.packages("umx")
# load as usual
library("umx")

# The current package version can be shown with:
umxVersion("umx")

# Get the latest NPSOL and multi-core build of OpenMx
install.OpenMx("NPSOL")

# Bleeding edge version of OpenMx for MacOS
install.OpenMx("travis")

# ============
# = CFA Code =
# ============

# Load the umx library (this is assumed in subsequent examples
library("umx")

# Load demo data consisting of 5 correlated variables, x1:x5 
data(demoOneFactor)

# Create a list of the manifest variables for use in specifying the model
manifests = paste0("x", 1:5) # 'x1', 'x2', ...'x5'

# Create model cfa1, with name 'CFA', data demoOneFactor, and the CFA paths.

cfa1 = umxRAM("CFA", data = demoOneFactor,
	# Create latent variable 'G', with fixed variance of 1 and mean of 0
	umxPath(v1m0 = "G"),
	# Create 5 manifest variables, x1:x5, with free variance and mean
	umxPath(v.m. = manifests),
	# Create 1-headed paths from G to each of the manifests
	umxPath("G", to = manifests)
)

# ====================
# = Parameter labels =
# ====================

x = xmuLabel(mxMatrix(name="means", "Full", ncol = 2, nrow = 2))
x$labels

# ========
# = Plot =
# ========

plot(cfa1, means = FALSE, fixed = TRUE)
plot(cfa1, std = TRUE, digits = 3, resid= 'line')


m1 = umxRAM("play", data = c("A", "B", "C"),
	umxPath(unique.pairs = c("A", "B", "C"))
)

# ==============================================
# = Inspecting model parameters and residuals. =
# ==============================================

# Show parameters, below .1, with label containing `x2' 
parameters(cfa1, "above", .5, pattern= "x2")

residuals(cfa1, suppress = .005)

# ==================================
# = Modifying and comparing models =
# ==================================

# Variable names in the Duncan data
dimnames = c("RespOccAsp", "RespEduAsp", "RespParAsp", "RespIQ", "RespSES",
             "FrndOccAsp", "FrndEduAsp", "FrndParAsp", "FrndIQ", "FrndSES")
# lower-triangle of correlations among these variables
tmp = c(
	0.6247,
	0.2137, 0.2742,
	0.4105, 0.4043, 0.1839,
	0.3240, 0.4047, 0.0489, 0.2220,
	0.3269, 0.3669, 0.1124, 0.2903, 0.3054,
	0.4216, 0.3275, 0.0839, 0.2598, 0.2786, 0.6404,
	0.0760, 0.0702, 0.1147, 0.1021, 0.0931, 0.2784, 0.1988, 
	0.2995, 0.2863, 0.0782, 0.3355, 0.2302, 0.5191, 0.5007,  0.2087,
	0.2930, 0.2407, 0.0186, 0.1861, 0.2707, 0.4105, 0.3607, -0.0438, 0.2950
)

# Use the umx_lower2full function to create a full correlation matrix
duncanCov = umx_lower2full(tmp, diag = FALSE, dimnames = dimnames)

# Turn the duncan data into an mxData object for the model
duncanCov = mxData(duncanCov, type = "cov", numObs = 300)

respondentFormants   = c("RespSES", "FrndSES", "RespIQ", "RespParAsp")
friendFormants       = c("FrndSES", "RespSES", "FrndIQ", "FrndParAsp")
latentAspiration     = c("RespLatentAsp", "FrndLatentAsp")
respondentOutcomeAsp = c("RespOccAsp", "RespEduAsp")
friendOutcomeAsp     = c("FrndOccAsp", "FrndEduAsp")


duncan1 = umxRAM("Duncan", data = duncanCov,
	# Working from the left of the model, as laid out in the figure, to right...

	# 1. Add all distinct paths between variables to allow the 
	# exogenous manifests to covary with each other.
	umxPath(unique.bivariate = c(friendFormants, respondentFormants)),

	# 2. Add variances for the exogenous manifests,
	# These are assumed to be error-free in this model,
	# and are fixed at their known value).
	umxPath(var = c(friendFormants, respondentFormants), fixedAt = 1),

	# 3. Paths from IQ, SES, and parental aspiration 
	# to latent aspiration for Respondents:
	umxPath(respondentFormants, to = "RespLatentAsp"),
	# And same for friends
	umxPath(friendFormants,     to = "FrndLatentAsp"),

	# 4. Add residual variance for the two aspiration latent traits.
	umxPath(var = latentAspiration),

	# 5. Allow the latent traits each influence the other.
	# This is done using fromEach, and the values are 
	# bounded to improve stability.
	# note: Using one-label would equate these 2 influences
	umxPath(fromEach = latentAspiration, lbound = 0, ubound = 1), 

	# 6. Allow latent aspiration to affect respondent's
	# occupational & educational aspiration.
	# note: firstAt = 1 is used to provide scale to the latent variables.
	umxPath("RespLatentAsp", to = respondentOutcomeAsp, firstAt = 1),

	# And their friends
	umxPath("FrndLatentAsp", to = friendOutcomeAsp, firstAt = 1),
	
	# 7. Finally, on the right hand side of figure, we add
	# residual variance for the endogenous manifests.
	umxPath(var = c(respondentOutcomeAsp, friendOutcomeAsp))
)

# ====================
# = Modifying models =
# ====================

# Collect a list of paths to drop
pathList = c("RespLatentAsp_to_FrndLatentAsp", "FrndLatentAsp_to_RespLatentAsp")

# Modify the model duncan1, requesting a comparison table:
duncan2 = umxModify(duncan1, update = pathList, name = "No_influence", comparison = TRUE)

# An example using regex, to drop all paths beginning "G_to_"
cfa2 = umxModify(cfa1, regex = "^G_to.*")

# ====================
# = Comparing models =
# ====================

umxCompare(duncan1, duncan2, report = "inline")

# To open the output as an html table in a browser, say:
umxCompare(duncan1, duncan2, report = "html")

# =============================
# = Equating model parameters =
# =============================

parameters(duncan1, pattern = "IQ_to_")


duncan3 = umxModify(duncan1, name = "Equate IQ effect", comparison = TRUE,
	master = "RespIQ_to_RespLatentAsp", 
	update = "FrndIQ_to_FrndLatentAsp" 
)

# ================
# = ACE examples =
# ================

require(umx); 
# open the built in dataset of Australian height and weight twin data
data("twinData")
selDVs = c("wt")
dz = twinData[twinData$zygosity == "DZFF", ]
mz = twinData[twinData$zygosity == "MZFF", ]

ACE1 = umxACE(selDVs = selDVs, dzData = dz, mzData = mz, sep = "")
ACE2 = umxModify(ACE1, update = "c_r1c1", name = "dropC")
umxSummary(ACE1, std = FALSE, report = 'html', digits = 3, comparison = ACE2)
parameters(ACE1)

ACE2 = umxModify(ACE1, update = "c_r1c1", name = "dropC")

# ================================
# = Example Common Pathway model =
# ================================

# load twin data built into umx
data("twinData")

# Selecting the 'ht' and 'wt' variables
selDVs = c("ht", "wt")
mzData = subset(twinData, zygosity == "MZFF",)
dzData = subset(twinData, zygosity == "DZFF",)

# Run and report a common-pathway model
CP1 = umxCP(selDVs = selDVs, dzData = dzData, mzData = mzData, suffix = "")

paths = c("c_cp_r1c1", "cs_r1c1", "cs_r2c2")
CP2 = umxModify(CP1, update = paths, name = "dropC", comparison = TRUE)

CP2 = umxModify(CP1, regex = "(^cs_)|(^c_cp_)", name = "dropC")
umxSummary(CP2, comparison = CP1)

# ====================================
# = Example Gene x environment model =
# ====================================

data("twinData")
twinData$age1 = twinData$age2 = twinData$age
# Define the DV and definition variables
selDVs  = c("bmi1", "bmi2")
selDefs = c("age1", "age2")
selVars = c(selDVs, selDefs)

# Create datasets
mzData  = subset(twinData, zygosity == "MZFF")
dzData  = subset(twinData, zygosity == "DZFF")

# Build, run and report the GxE model using selected DV and moderator
# umxGxE will remove and report rows with missing data in definition variables.
GE1 = umxGxE(selDVs = selDVs, selDefs = selDefs, 
  dzData = dzData, mzData = mzData,  dropMissingDef = TRUE)

# Shift the legend to the top right
umxSummary(GE1, location = "topright")

# plot standardized and raw output in separate graphs
umxSummary(GE1, separateGraphs = TRUE)

GE2 = umxModify(GE1, update = "am_r1c1", comparison = TRUE)
umxReduce(GE1)

# =================================
# = Example GxE windowed analysis =
# =================================

require(umx);
data("twinData") 
mod     = "age"
selDVs  = c("bmi1", "bmi2")

# select the younger cohort of twins
tmpTwin = twinData[twinData$cohort == "younger", ]
# Drop twins with missing moderator
tmpTwin = tmpTwin[!is.na(tmpTwin[mod]), ]
mzData  = subset(tmpTwin, zygosity == "MZFF", c(selDVs, mod))
dzData  = subset(tmpTwin, zygosity == "DZFF", c(selDVs, mod))
# toggle autoplot off, so we don't plot every level of the moderator
umx_set_auto_plot(FALSE)
umxGxE_window(selDVs = selDVs, moderator = mod, mzData = mzData, dzData = dzData)
umx_set_auto_plot(TRUE)


## End(Not run)

Get the expected vcov matrix

Description

Extract the expected covariance matrix from an OpenMx::mxModel()

Usage

umxExpCov(object, latents = FALSE, manifests = TRUE, digits = NULL, ...)

Arguments

Value

• expected covariance matrix

See Also

• umxRun(), umxCI_boot()

Other Reporting functions: RMSEA(), RMSEA.MxModel(), RMSEA.summary.mxmodel(), extractAIC.MxModel(), loadings(), loadings.MxModel(), residuals.MxModel(), tmx_show(), tmx_show.MxMatrix(), umxCI(), umxCI_boot(), umxConfint(), umxExpMeans(), umxFitIndices(), umxRotate()

Examples

## Not run: 
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)
)#'
vcov(m1) # supplied by OpenMx
umxExpCov(m1, digits = 3)


## End(Not run)

Extract the expected means matrix from an OpenMx::mxModel()

Description

Extract the expected means matrix from an OpenMx::mxModel()

Usage

umxExpMeans(model, manifests = TRUE, latents = NULL, digits = NULL)

Arguments

Value

• expected means

See Also

Other Reporting functions: RMSEA(), RMSEA.MxModel(), RMSEA.summary.mxmodel(), extractAIC.MxModel(), loadings(), loadings.MxModel(), residuals.MxModel(), tmx_show(), tmx_show.MxMatrix(), umxCI(), umxCI_boot(), umxConfint(), umxExpCov(), umxFitIndices(), umxRotate()

Examples

## Not run: 
require(umx)
data(demoOneFactor)
manifests = names(demoOneFactor)
m1 = umxRAM("One Factor", data = demoOneFactor,
	umxPath("G", to = manifests),
	umxPath(v.m. = manifests),
	umxPath(v1m0 = "G")
)

umxExpMeans(m1)
umxExpMeans(m1, digits = 3)

## End(Not run)

umxFactor

Description

A convenient version of OpenMx::mxFactor() supporting the common case in which the factor levels are those in the variable.

Usage

umxFactor(
  x = character(),
  levels = NULL,
  labels = levels,
  exclude = NA,
  ordered = TRUE,
  collapse = FALSE,
  verbose = FALSE,
  sep = NA
)

Arguments

Value

• OpenMx::mxFactor()

References

• https://github.com/tbates/umx, https://tbates.github.io

See Also

• umxFactanal(), OpenMx::mxFactor()

Other Data Functions: noNAs(), prolific_anonymize(), prolific_check_ID(), prolific_read_demog(), umx, umxHetCor(), umx_as_numeric(), umx_cont_2_quantiles(), umx_lower2full(), umx_make_MR_data(), umx_make_TwinData(), umx_make_fake_data(), umx_make_raw_from_cov(), umx_merge_randomized_columns(), umx_polychoric(), umx_polypairwise(), umx_polytriowise(), umx_read_lower(), umx_rename(), umx_reorder(), umx_score_scale(), umx_select_valid(), umx_stack(), umx_strings2numeric()

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)

Return factor scores from a model as an easily consumable dataframe.

Description

umxFactorScores takes a model, and computes factors scores using the selected method (one of 'ML', 'WeightedML', or 'Regression') It is a simple wrapper around mxFactorScores. For missing data, you must specify the least number of variables allowed for a score (subjects with fewer than minManifests will return a score of NA.

Usage

umxFactorScores(
  model,
  type = c("ML", "WeightedML", "Regression"),
  minManifests = NA,
  return = c("Scores", "StandardErrors")
)

Arguments

Value

• dataframe of scores.

References

• https://github.com/tbates/umx, https://tbates.github.io

See Also

• OpenMx::mxFactorScores()

Other Reporting Functions: umx, umxAPA(), umxGetLatents(), umxGetManifests(), umxGetModel(), umxGetParameters(), umxParameters(), umx_aggregate(), umx_time()

Examples

## Not run: 
m1 = umxEFA(mtcars, factors = 2)
x = umxFactorScores(m1, type = 'Regression', minManifests = 3)

# =========================================================================
# = histogram of F1 and plot of F1 against F2 showing they are orthogonal =
# =========================================================================
hist(x$F1)
plot(F1 ~ F2, data = x)

m1 = umxEFA(mtcars, factors = 1)
x = umxFactorScores(m1, type = 'Regression', minManifests = 3)
x

## End(Not run)

Get additional fit-indices for a model with umxFitIndices

Description

Computes a variety of fit indices.

Usage

umxFitIndices(model, ...)

Arguments