Multilevel Models for Twins Nested in Pairs Using nlme::lme

Separate means and variance components by zygosity · Last updated: October 2, 2026
R script: Twin_nlme_Example.R (uses EDF9850_Functions.R)

1. Design and data

Twins (level 1) are nested within pairs (level 2). Monozygotic (MZ) and dizygotic (DZ) pairs are expected to differ in how similar the twins are, so the means, the pair-level (between-pair) variance, and the twin-level (within-pair, residual) variance are each estimated separately by zygosity. The example data are simulated (300 MZ and 300 DZ pairs; 1200 twins) with known variances: pair variances of 6 (MZ) and 3 (DZ) and residual variances of 2 (MZ) and 5 (DZ), before the effects of smoking and age, which are part of the variance in the empty model.

• Outcome: y (symptoms)
• Twin predictor: smoke (0 = no, 1 = yes), which varies within pairs
• Pair predictor: age, which is identical for both twins in a pair
• IDs: pairID (the level-2 cluster) and twinID (1 or 2 within each pair)

2. Predictor setup

Zygosity binary contrasts are needed to give each group its own random-intercept variance. Age is centered at 40. Smoking varies both within and between pairs, so it is split into a between-pair part and a within-pair part, so that each effect is interpreted at its own level: addUnitMeans() (from EDF9850_Functions.R) adds the pair (cluster) mean CMsmoke, which is centered at .40 (close to its overall mean of .403) as CMsmoke40, and WCsmoke is each twin’s deviation from the pair mean.

twins$MZ = as.numeric(twins$zyg == "MZ")
twins$DZ = as.numeric(twins$zyg == "DZ")
twins$age40 = twins$age - 40
twins = addUnitMeans(data=twins, unitVariable="pairID",
                     meanVariables=c("smoke"), newNames=c("CMsmoke"))
twins$CMsmoke40 = twins$CMsmoke - .40       # between-pair (pair mean)
twins$WCsmoke   = twins$smoke - twins$CMsmoke   # within-pair (twin deviation)

3. Key syntax for zygosity-specific variances

random=list(pairID=pdBlocked(list(pdIdent(~0+MZ),pdIdent(~0+DZ))))
weights = varIdent(form=~1|zyg)

• Pair level: pdBlocked with the two binary contrasts creates two random intercepts, one per zygosity, each with its own variance. Each is zero for the other group’s pairs. A pair variance is the covariance between co-twins.
• Twin level: varIdent lets the residual SD differ by zygosity. The MZ SD is the reference (sigma); the DZ SD is a multiple of it.
• lmeSummary() (from EDF9850_Functions.R) prints these models directly: the pair variances (the MZ and DZ random effects), a “Residual Variance by Group” table (the plain “Residual” row is only the MZ reference value), and, on request, ICC1, ICC2, design effect, and random-intercept CIs for each zygosity.
• All models use method="REML", including the likelihood ratio test of variance structures (the fixed effects are identical across those models). Sets of fixed effects are tested jointly with conMDSummary() (section 6). The control=ctrl argument (lmeControl(maxIter=200, msMaxIter=200, opt="optim")) only makes convergence more stable; these models also converge with the defaults.

4. Model 1: empty model with means and variances by zygosity

y ~ 0 + zyg removes the common intercept so the two fixed effects are the MZ and DZ means directly (the same “means per zygosity” parameterization as NOINT in SAS).

Model1 = lme(data=twins, method="REML", control=ctrl,
         fixed=y~0+zyg,
         random=list(pairID=pdBlocked(list(pdIdent(~0+MZ),pdIdent(~0+DZ)))),
         weights = varIdent(form=~1|zyg))
s1=lmeSummary(Model1, CI=TRUE, ICC1=TRUE, ICC2=TRUE, randCI=TRUE, LRT=TRUE)

This one call prints the means by zygosity (with CIs), the pair and residual variances by zygosity, ICC1 and ICC2 by zygosity, the random-intercept CIs by zygosity, and an LRT against a model with no random effects (REML, with the same residual variance function; χ²(2) = 302.5, p < .001). The saved result s1 is reused below (s1$ICC1 holds the pair variance, residual variance, and ICC1 for each zygosity).

ICC1 is the pair variance divided by total variance within each zygosity (the twin correlation); ICC2 is the reliability of the pair means, which for two twins per pair uses the residual variance divided by 2.

Do the variances differ by zygosity? Compare Model 1 to a model with one pair variance and one residual variance for both zygosities. That model, Model1eq, has the same fixed effects but a single random intercept (~ 1 | pairID) and no weights argument, so it estimates only one pair variance and one residual variance. Because the fixed effects are identical in the two models, the REML likelihood ratio test is valid and no ML refit is needed:

Model1eq = lme(data=twins, method="REML", control=ctrl,
               y~0+zyg, random=~1 | pairID)
randLRT(Model1eq, Model1)
    Model Npar    LL -2LL  AIC  BIC    LRT DF p.value p.mix
 Model1eq    4 -2978 5956 5964 5984
   Model1    6 -2961 5921 5933 5964 34.743  2   <.001 <.001

Model1eq has 4 parameters (two means, one pair variance, one residual variance) and Model1 has 6 (two means, two pair variances, two residual variances), so the test has 2 df. randLRT() (from EDF9850_Functions.R) compares the two REML fits directly, taking the df from the difference in parameters. Use its p.value column here: p.mix is the boundary-corrected p-value for testing whether a variance is zero, but these variances are tested for being equal, not zero, so plain χ²(2) is the right reference distribution. Allowing separate variances fits much better (χ²(2) = 34.7, p < .001).

5. Model 2: add smoking (twin) and age (pair)

Model2 = lme(data=twins, method="REML", control=ctrl,
         fixed=y~0+zyg+WCsmoke+CMsmoke40+age40,
         random=list(pairID=pdBlocked(list(pdIdent(~0+MZ),pdIdent(~0+DZ)))),
         weights = varIdent(form=~1|zyg))
s2=lmeSummary(Model2)

lmeSummary() prints the fixed effects, the conditional variances by zygosity, and an overall test that all slopes are zero. That overall test is an F-test with Satterthwaite denominator df (the slopes are at different levels, twin-level smoking and pair-level age, so nlme’s own anova(Terms = ) could not combine them): F(3, 541.4) = 129.03, p < .001. The two zygosity means are left out of this test, since testing that they equal zero would be meaningless.

The proportional reduction in each variance component relative to Model 1 (a pseudo-R²) comes from PseudoR2(). For an lme model with varIdent it gives one residual row per zygosity, each compared to the same group in Model 1:

PseudoR2(Model1, Model2)
 Term           Model  Variance Pseudo.R2
 MZ             Model1    7.857
 MZ             Model2    6.937     0.117
 DZ             Model1    4.710
 DZ             Model2    2.745     0.417
 Residual (MZ)  Model1    2.868
 Residual (MZ)  Model2    1.914     0.332
 Residual (DZ)  Model1    5.492
 Residual (DZ)  Model2    4.764     0.133

Interpretation. The within-pair smoking effect compares a smoking twin to the co-twin average in the same pair; smokers score 2.03 higher, controlling for everything shared within the pair. The between-pair effect (2.30) compares pairs that differ by one unit in the proportion of twins who smoke. Age is a pair predictor, so it can reduce only pair-level variance; smoking has a within-pair part, so it can reduce twin-level variance too.

The means (intercepts) are separate by zygosity, as in Model 1; the slopes are shared across zygosities.

6. Joint tests of sets of fixed effects

All models are estimated with REML; no ML refits are used. Sets of fixed effects are tested jointly with conMDSummary() (from EDF9850_Functions.R), a test of L b = 0, where b is the vector of fixed effects and each row of L is one linear combination of them. For an lme model it is an F-test with Satterthwaite denominator df: the same df method lmeSummary() uses for its Fixed Effects table, extended to several df at once the way lmerTest::contestMD() does for lmer. (nlme itself offers only “containment” df, and anova(Terms = ) cannot combine terms with different denominator df, such as twin-level smoking and pair-level age.) The Satterthwaite df for lme are approximate, so they can differ slightly from the same model in lmer, and they are computed with a fixed random seed so they are reproducible. With ddf = "asymptotic" the test is instead a multivariate Wald χ², W = (Lb)′ [L V L′]⁻¹ (Lb), with df equal to the number of rows of L, which needs no denominator df.

• conMDSummary(model, terms=c(...)) tests that a subset of slopes are all zero.
• conMDSummary(model, L=...) tests any linear combinations. L can be a character vector of fixed-effect names (each tested = 0), a numeric vector named after the fixed effects (unnamed effects are 0, e.g. c(zygDZ=1, zygMZ=-1) for DZ minus MZ), or a list mixing those forms (each numeric vector is one row, labeled by its list name if it has one; each character element adds one row per fixed-effect name, so a difference can simply be followed by plain coefficient names). It works for lme, gls, lmer and glmer models.

To ask whether the slopes differ by zygosity, Model3 adds the interactions of zyg with each predictor to Model 2. The means stay separate by zygosity, and MZ is the reference group for the interactions, so zygDZ:WCsmoke is the DZ-minus-MZ difference in the WCsmoke slope.

Model3 = lme(data=twins, method="REML", control=ctrl,
         fixed=y~0+zyg+WCsmoke+CMsmoke40+age40+zyg:WCsmoke+zyg:CMsmoke40+zyg:age40,
         random=list(pairID=pdBlocked(list(pdIdent(~0+MZ),pdIdent(~0+DZ)))),
         weights = varIdent(form=~1|zyg))
conMDSummary(Model2, terms=c("WCsmoke","CMsmoke40","age40"))
conMDSummary(Model3, L=c("zygDZ:WCsmoke","zygDZ:CMsmoke40","zygDZ:age40"))
conMDSummary(Model3, L=list("DZ - MZ mean"=c(zygDZ=1, zygMZ=-1),
                            c("zygDZ:WCsmoke","zygDZ:CMsmoke40","zygDZ:age40")))
conMDSummary(Model2, L=c(WCsmoke=1, CMsmoke40=-1))
conMDSummary(Model2, L=c(zygDZ=1, zygMZ=-1))

The interactions are not significant as a set, so the slopes can be treated as equal across zygosities, and the means do not differ either. The within- and between-pair smoking effects (2.03 and 2.30) also do not differ significantly.

7. ACE point estimates

The pair-level variance in each zygosity is the covariance between co-twins, so ACE components can be derived from Model 1 (the empty model; the decomposition applies to the total variance):

cov(MZ) = A + C        cov(DZ) = ½A + C        residual variance (MZ) = E
A = 2 (covMZ − covDZ)    C = 2 covDZ − covMZ    E = residual variance

a2 = 2 (rMZ − rDZ)    c2 = 2 rDZ − rMZ    e2 = 1 − rMZ      (r = ICC1; Falconer)

This model does not constrain total variance to be equal in MZ and DZ twins, and the estimates are not bounded at zero (A or C can be negative in a sample).

aceFromTable = function(icc, digits=3) {   # icc = s1$ICC1 from lmeSummary(..., ICC1=TRUE)
  covMZ = icc$Cluster.Var[1]; covDZ = icc$Cluster.Var[2]
  rMZ   = icc$ICC1[1];        rDZ   = icc$ICC1[2]
  round(c(A = 2 * (covMZ - covDZ), C = 2 * covDZ - covMZ, E = icc$Resid.Var[1],
          a2 = 2 * (rMZ - rDZ),    c2 = 2 * rDZ - rMZ,    e2 = 1 - rMZ), digits)
}
aceFromTable(s1$ICC1)
    A     C     E    a2    c2    e2
6.294 1.563 2.868 0.542 0.191 0.267

The estimates are rounded to digits decimal places (default 3).

8. Delta-method standard errors and CIs

lme does not estimate A, C, E, a2, c2 and e2 directly, but it stores the covariance matrix of its variance parameters in m$apVar. The parameters are on the log scale: log pair SD (MZ), log pair SD (DZ), log DZ/MZ residual SD ratio, and log residual SD (MZ). The delta method gives Var(g(θ)) ≈ J V J′, where V is that matrix and J contains the derivatives of each quantity with respect to the parameters. This is what the SAS NLMIXED ESTIMATE statements do. The script computes J numerically, so no extra package is needed, and checks that the function reproduces the direct point estimates. aceDelta(m, level=0.95, digits=3) rounds its table to digits decimal places.

aceDelta(Model1)

The SEs and CIs for the MZ and DZ means are in the lmeSummary(Model1, CI=TRUE) output. The ACE intervals here are symmetric Wald intervals (normal critical value; NLMIXED uses t). They can fall outside the possible range of a variance or proportion (as the interval for C does here) and are only approximate for variance-type quantities, especially A and C. Heritability is substantial (a2 ≈ .54), and the shared-environment component is small, with an interval for c2 that only just excludes zero.