Fits classical biometric variance–decomposition models using a linear mixed model formulation implemented with nlme::lme.
Usage
ACDE(model, data, type = c("AE", "ACE", "ADE"), method = "ML")Value
Object of class "ACDEfit" containing:
fit
nlme::lmeobjectvar Variance components (A,C,D,E)
h2 Narrow-sense heritability
c2 Shared environment (ACE only)
d2 Dominance (ADE only)
H2 Broad-sense heritability (ADE only)
Details
The function estimates additive genetic (A), shared environmental (C), dominance genetic (D), and unique environmental (E) variance components from nuclear family data (parents and offspring) in long format.
Supported family structures:
Parent-child trios (one offspring)
Nuclear families with any number of siblings
The function automatically detects the family structure and selects the correct additive-genetic parameterisation.
Only AE, ACE and ADE models are fitted.
Phenotypic variance is decomposed as $$V_P = V_A + V_C + V_D + V_E$$
The model is fitted as a linear mixed model with family-level random effects and individual residual variance.
Automatic family detection
The function detects whether families contain one or multiple offspring.
Trios: collapsed additive parameterisation.
Siblings: transmission decomposition parameterisation.
Trio parameterisation
Additive effects represented as: $$A = 0.5 Mother + 0.5 Father + 1 Child$$
This reproduces the expected parent–offspring covariance: $$Cov = 1/2 V_A$$
Multi-sibling parameterisation
Additive genetic variance is decomposed into:
maternal transmission (\(A_m\))
paternal transmission (\(A_f\))
Mendelian sampling (\(M_s\))
For offspring: $$A = A_m + A_f + M_s$$
Total additive variance: $$V_A = 2(\sigma^2_{Am} + \sigma^2_{Af}) + \sigma^2_{Ms}$$
This produces correct covariances:
Parent–offspring: \(1/2 V_A\)
Sibling–sibling: \(1/2 V_A\)
This formulation generalises to any number of siblings.
Required columns in data
familyid Nuclear family identifier.
var1 Maternal transmission coefficient.
var2 Paternal transmission coefficient.
var3 Offspring (Mendelian sampling) indicator.
These variables encode expected genetic transmission and are not role indicators.
Coding for a nuclear family:
role var1 var2 var3
-------------------------
father 0 1 0
mother 1 0 0
child 1 1 1
For multiple siblings, each offspring receives identical coding:
role var1 var2 var3
-------------------------
father 0 1 0
mother 1 0 0
sib1 1 1 1
sib2 1 1 1
sib3 1 1 1
Examples
library(nlme)
set.seed(1)
simulate_families <- function(n_fam = 200)
{
VA <- 0.4; VC <- 0.2; VD <- 0.1; VE <- 0.3
out <- list()
for(f in 1:n_fam){
Cfam <- rnorm(1,0,sqrt(VC))
Af <- rnorm(1,0,sqrt(VA))
Am <- rnorm(1,0,sqrt(VA))
make_child <- function(){
Mend <- rnorm(1,0,sqrt(0.5*VA))
A <- 0.5*(Af+Am)+Mend
D <- rnorm(1,0,sqrt(VD))
E <- rnorm(1,0,sqrt(VE))
A + D + Cfam + E
}
out[[f]] <- data.frame(
familyid=f,
role=c("father","mother","sib1","sib2","sib3"),
y=c(
Af + Cfam + rnorm(1,0,sqrt(VE)),
Am + Cfam + rnorm(1,0,sqrt(VE)),
make_child(), make_child(), make_child()
)
)
}
dat <- do.call(rbind,out)
dat$var1 <- as.integer(dat$role=="mother")
dat$var2 <- as.integer(dat$role=="father")
dat$var3 <- as.integer(grepl("sib", dat$role))
dat
}
dat <- simulate_families()
AE <- ACDE(y~1, dat, "AE")
ACE <- ACDE(y~1, dat, "ACE")
ADE <- ACDE(y~1, dat, "ADE")
anova(AE$fit, ACE$fit, ADE$fit)
#> Model df AIC BIC logLik Test L.Ratio p-value
#> AE$fit 1 5 2772.272 2796.811 -1381.136
#> ACE$fit 2 6 2637.652 2667.099 -1312.826 1 vs 2 136.6198 <.0001
#> ADE$fit 3 6 2637.652 2667.099 -1312.826
############################################################
# Create rectangular variance table (important!)
############################################################
summary(AE)
#>
#> Summary of AE model
#> -----------------------------------
#> Parameterisation: sibling
#>
#> Variance components:
#> A C D E
#> 1.493 0.000 0.000 0.656
#>
#> Variance proportions:
#> h2
#> 0.695
#>
#> Total phenotypic variance: 2.149
#> logLik: -1381.136
summary(ACE)
#>
#> Summary of ACE model
#> -----------------------------------
#> Parameterisation: sibling
#>
#> Variance components:
#> A C D E
#> 0.000 0.399 0.000 0.604
#>
#> Variance proportions:
#> h2 c2
#> 0.000 0.398
#>
#> Total phenotypic variance: 1.003
#> logLik: -1312.826
summary(ADE)
#>
#> Summary of ADE model
#> -----------------------------------
#> Parameterisation: sibling
#>
#> Variance components:
#> A C D E
#> 0.000 0.000 0.399 0.604
#>
#> Variance proportions:
#> h2 d2 H2
#> 0.000 0.398 0.398
#>
#> Total phenotypic variance: 1.003
#> logLik: -1312.826
require(gap.datasets)
model <- bwt ~ male + first + midage + highage + birthyr
AE <- ACDE(model,mfblong)
ACE <- ACDE(model,mfblong,type="ACE")
ADE <- ACDE(model,mfblong,type="ADE")
anova(AE$fit,ACE$fit,ADE$fit)
#> Model df AIC BIC logLik Test L.Ratio p-value
#> AE$fit 1 8 45539.51 45587.56 -22761.76
#> ACE$fit 2 9 45535.98 45590.04 -22758.99 1 vs 2 5.530908 0.0187
#> ADE$fit 3 9 45535.98 45590.04 -22758.99