umxACE_GLM: Cholesky ACE twin model with mxFitFunctionGLM
In umx: Structural Equation Modeling and Twin Modeling in R

umxACE_GLM

R Documentation

Cholesky ACE twin model with mxFitFunctionGLM

Description

RAM MZ/DZ Cholesky ACE for exponential-family phenotypes. Not a patch on umxACE(): that function uses mxExpectationNormal and a %*% t(a) algebras, which GLM cannot use.

Additive and shared factors have variance 1. MZ uses one shared a_i and (if dzCr is 1) one shared c_i per factor: two copies with correlation 1 are singular for Laplace. DZ has a_i_T1 with a_i_T2 covariance dzAr. Unique e is always per twin. Loadings use umxACE labels ⁠a_r{row}c{col}⁠. Off-diagonal e is the GLM analogue of E covariances (shared unique / measurement error). Family leftover is still independent given eta. Unique environment is the e latent and the family leftover. nAGQ is 1. Cholesky is nVar-general. MZ Laplace dimension is 4 * nVar latents.

Bernoulli: fix e diagonal at 1. Leftover (probit 1, logit \pi^2/3) already sets the observation scale. Univariate ACE still has three variance components and only two twin correlations, so a free e_ii is not identified (e went to 0 and a was eaten). For binomial items, e_rjc_j is fixed at 1 (same idea as fixing unique scale in threshold ACE). Off-diagonal e stays free so multivariate unique covariance is not forced to 0. Do not set e11 to leftover: that would double-count unique on eta. umxSummaryACE_GLM() still uses ⁠E = ee' + D⁠ for shares. Poisson keeps a free e diagonal (leftover is on the counts). Poisson ACE starts the Cholesky diagonals at a=\sqrt{0.8}, e=\sqrt{0.2}, c at 0.2 (80/20 A/E on eta; C not pinned to the floor) so the optimizer walks downhill from a genetic model. Free ACE diagonals have lbound = 1e-4 so the derivative stays alive at the floor.

Gamma and inverse Gaussian: not for ACE. Inverse Gaussian is the hitting-time law of Brownian motion with drift (reaction time, age of onset). Gamma is the constant-CV amount / duration model. Both need trial-level data and, for a process interpretation, structure on two parameters (mean and shape/lambda, or drift and barrier with diffusion scale fixed). Twin ACE is one phenotype per person, so those uses are out of scope. The wrapper still accepts the families: leftover shape (Gamma) or lambda (inverse Gaussian) is fixed at 1 (phi = 1/theta) unless you pass a numeric theta; a free theta label is rejected (collides with e). That is not the cognitive-psych constraint of fixing diffusion scale \sigma. Leave them in the too-hard box; use Poisson for counts and threshold umxACE() for binary.

Binomial ⁠trials=⁠ is not implemented. mxFamily() can take trials as a raw-data column (n successes out of n_i trials per row). This wrapper does not: there is no trials argument, and MZ/DZ data are subset to the phenotype ⁠_T1⁠/⁠_T2⁠ columns, so any n columns are dropped. Binomial here is Bernoulli (n = 1) only. Grouped or item-sum binomials with a dynamic n cannot be wired. Most real binomial counts are overdispersed; there is no beta-binomial (or other extra-binomial leftover) to absorb that, so users will keep using liability umxACE() with more than one threshold, or a Gaussian for large counts. Negative binomial on the counts is the overdispersed alternative when there is no upper bound.

Hessian / standard errors stay off. umxRAM_GLM() (no latents) yields finite SEs. Twin ACE always has several latents. On univariate Poisson, n = 50 pairs gave a Hessian that was not positive definite and SEs all NA (optimizer status 6). At n = 200 the Hessian was positive definite and SEs were finite, but SEs on a and c were about 1 on loadings of about 0.5, and status was still 6. Default remains ⁠Calculate Hessian⁠ / ⁠Standard Errors⁠ = No. After a fitted model, family-score SEs: mxComputeOnce('fitfunction', 'information', 'meat') then mxComputeStandardError(). Poisson intercept matches glm(); Poisson ACE intercept matches a pair bootstrap in tests.

Usage

umxACE_GLM(
  name = "ACE_GLM",
  selDVs,
  mzData,
  dzData,
  sep = "_T",
  family = stats::poisson(),
  theta = NULL,
  nAGQ = 1L,
  dzAr = 0.5,
  dzCr = 1,
  equateMeans = TRUE,
  autoRun = getOption("umx_auto_run"),
  tryHard = c("no", "yes", "ordinal", "search")
)

Arguments

`name`	model name (default "ACE_GLM")
`selDVs`	base phenotype names (sep expands to `⁠_T1⁠`/`⁠_T2⁠`). Cholesky is `nVar`-general.
`mzData`	MZ data.frame
`dzData`	DZ data.frame
`sep`	twin separator (default `"_T"`; must not be `""`)
`family`	`stats::family` or `mxFamily()` for the phenotype (both twins). Poisson is the intended ACE family. Binomial is Bernoulli only: `⁠trials=⁠` and per-row `n` columns are not implemented (see Details). Scale-identified (`e` diagonal fixed at 1) but leftover eats heritability — use threshold `umxACE()` for binary. Gamma and inverse Gaussian are accepted but not supported for ACE (trial-level process models; see Details).
`theta`	passed to `mxFamily()`. For Gamma, leftover shape; for inverse Gaussian, leftover lambda. Both fixed (default 1, `phi = 1`). A free-parameter label is rejected (collides with `e`). Not used for Poisson or binomial.
`nAGQ`	must be 1 (more than one latent)
`dzAr`	DZ additive genetic correlation (default 0.5)
`dzCr`	DZ C correlation (default 1; 0.25 for ADE)
`equateMeans`	equate T1/T2 intercepts (default TRUE)
`autoRun`	whether to run
`tryHard`	`xmu_safe_run_summary()` tryHard

Value

supermodel of class MxModelACE_GLM

Examples

## Not run: 
set.seed(1)
n = 40
a = 0.4; c = 0.2; e = 0.5; b0 = 0.3
simPair = function(rA, n) {
	A = rnorm(n)
	C = rnorm(n)
	E1 = rnorm(n)
	E2 = rnorm(n)
	A2 = rA * A + sqrt(1 - rA^2) * rnorm(n)
	y1 = rpois(n, lambda = exp(b0 + a * A + c * C + e * E1))
	y2 = rpois(n, lambda = exp(b0 + a * A2 + c * C + e * E2))
	data.frame(y_T1 = y1, y_T2 = y2)
}
mzData = simPair(1, n)
dzData = simPair(0.5, n)
m1 = umxACE_GLM(selDVs = "y", mzData = mzData, dzData = dzData, sep = "_T")

## End(Not run)

umx documentation built on Oct. 11, 2026, 5:07 p.m.