Analysing within- and between-subject effects with semnova

library(semnova)

This vignette walks through everything the package currently does, using the data that ships with it. Every feature that is implemented appears here; if a function is exported and meant for users, you will find it below.

What semnova is for

semnova fits repeated-measures ANOVA designs as structural equation models. The advantage over a classical RM-ANOVA is that the dependent variable, and any covariates, can be latent. They are measured by several indicators, with measurement error taken out, rather than being manifest composites.

Under the hood, each design cell becomes a latent variable (eta), those are transformed by a contrast matrix into contrast variables (pi), and every effect you test is a general linear hypothesis on the pi layer. The maths is in vignette("glh", package = "semnova"); you do not need it to use the package.

The example data

Two versions of a reading study ship with the package: reading_manifest (one row per person and design cell) and reading_latent (the long form, with three indicators per cell). The design is 3 (grade) x 2 (sentence type), within subjects, with sex available as a between-subjects factor and iq as a covariate.

str(reading_manifest[c("id", "grade", "sentence", "sex", "re_fix_time", "iq")])
#> Classes 'tbl_df', 'tbl' and 'data.frame':    1602 obs. of  6 variables:
#>  $ id         : Named int  1 1 1 1 1 1 2 2 2 2 ...
#>   ..- attr(*, "names")= chr [1:1602] "002SG21032005" "002SG21032005" "002SG21032005" "002SG21032005" ...
#>  $ grade      : Factor w/ 6 levels "grade0","grade1",..: 2 3 5 2 3 5 2 3 5 2 ...
#>  $ sentence   : Factor w/ 2 levels "sentence","landolt": 1 1 1 2 2 2 1 1 1 2 ...
#>  $ sex        : chr  "m" "m" "m" "m" ...
#>  $ re_fix_time: num  NA 1.44 0.405 0.196 0.251 ...
#>  $ iq         : num  111 111 111 111 111 ...
table(reading_latent$indicator)
#> 
#>    fix_count ini_land_pos  re_fix_time 
#>         1602         1602         1602

A manifest within-subjects design

The simplest case: a single observed dependent variable, two crossed within-subject factors. Give semnova() the long-format data, the case identifier, the dependent variable, and the within factors.

contrasts_arg controls how each factor is coded. Anything R’s own contrast functions provide will do. Here we use sum-to-zero contrasts for both factors.

fit <- semnova(
    data = reading_manifest,
    id = "id",
    dv = "re_fix_time",
    withins = c("grade", "sentence"),
    contrasts_arg = list(grade = "contr.sum", sentence = "contr.sum")
)

Printing the fit gives you a description of the model and the ANOVA table: one row per model term, tested as a general linear hypothesis.

fit
#> semnova fit
#>   dependent variable : manifest
#>   occasions          : 6
#>   fit                : saturated (df = 0)
#> 
#> Hypothesis tests:
#>                   Chisq Df Pr(>Chisq)    
#> (Intercept)      4.7107  1    0.02998 *  
#> grade          832.0658  2    < 2e-16 ***
#> sentence       515.3166  1    < 2e-16 ***
#> grade:sentence 234.9193  2    < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova() on a single fit returns that table on its own, so you can work with it programmatically. Given several fits it compares them as nested models instead, which the invariance and covariance-structure sections below use.

anova(fit)
#> Hypothesis tests:
#>                   Chisq Df Pr(>Chisq)    
#> (Intercept)      4.7107  1    0.02998 *  
#> grade          832.0658  2    < 2e-16 ***
#> sentence       515.3166  1    < 2e-16 ***
#> grade:sentence 234.9193  2    < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova_table() is the same thing under a name that says what it returns, for when you would rather not go through the generic.

identical(anova_table(fit), anova(fit))
#> [1] TRUE

summary() adds the estimates underneath each omnibus test.

summary(fit)
#> semnova fit
#>   dependent variable : manifest
#>   occasions          : 6
#>   fit                : saturated (df = 0)
#> 
#> Hypothesis tests:
#>                   Chisq Df Pr(>Chisq)    
#> (Intercept)      4.7107  1    0.02998 *  
#> grade          832.0658  2    < 2e-16 ***
#> sentence       515.3166  1    < 2e-16 ***
#> grade:sentence 234.9193  2    < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Contrast estimates:
#>                                  estimate std_error z value ci_lower ci_upper
#> (Intercept)                       -0.0677   0.03121  -2.170  -0.1289  -0.0066
#> grade: grade1                      0.5915   0.02330  25.384   0.5458   0.6372
#> grade: grade2                     -0.0557   0.01636  -3.404  -0.0877  -0.0236
#> sentence: sentence1                0.5444   0.02398  22.701   0.4974   0.5914
#> grade:sentence: grade1:sentence1   0.2781   0.01991  13.971   0.2391   0.3171
#> grade:sentence: grade2:sentence1  -0.0272   0.01441  -1.886  -0.0554   0.0011
#>                                  sd_between Pr(>|z|)    
#> (Intercept)                          0.4045 0.029976 *  
#> grade: grade1                        0.3020  < 2e-16 ***
#> grade: grade2                        0.2120 0.000664 ***
#> sentence: sentence1                  0.3109  < 2e-16 ***
#> grade:sentence: grade1:sentence1     0.2580  < 2e-16 ***
#> grade:sentence: grade2:sentence1     0.1867 0.059276 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Adding a between-subjects factor

Naming a betweens factor fits a multigroup model, one group per level. It extends the ANOVA table with the between-subject term and every interaction involving it.

fit_between <- semnova(
    data = reading_manifest,
    id = "id",
    dv = "re_fix_time",
    withins = c("grade", "sentence"),
    betweens = "sex",
    contrasts_arg = list(
        grade = "contr.sum",
        sentence = "contr.sum",
        sex = "contr.treatment"
    )
)

anova_table(fit_between)
#> Hypothesis tests:
#>                       Chisq Df Pr(>Chisq)    
#> (Intercept)          1.5049  1     0.2199    
#> grade              354.7554  2     <2e-16 ***
#> sentence           381.8822  1     <2e-16 ***
#> grade:sentence      84.3705  2     <2e-16 ***
#> sex                  0.1923  1     0.6610    
#> grade:sex            1.9687  2     0.3737    
#> sentence:sex         1.8005  1     0.1797    
#> grade:sentence:sex   1.9587  2     0.3756    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A latent dependent variable

To make the dependent variable latent, point indicator at the column that says which indicator each row holds. Every design cell is then measured by all of the indicators, and the measurement model is estimated along with everything else.

We use two of the three available indicators here.

reading_two <- reading_latent[reading_latent$indicator != "re_fix_time", ]

fit_latent <- semnova(
    data = reading_two,
    id = "id",
    dv = "dv",
    indicator = "indicator",
    withins = c("grade", "sentence"),
    contrasts_arg = list(grade = "contr.sum", sentence = "contr.sum")
)

anova_table(fit_latent)
#> Hypothesis tests:
#>                   Chisq Df Pr(>Chisq)    
#> (Intercept)      6.8854  1    0.00869 ** 
#> grade          792.7875  2    < 2e-16 ***
#> sentence       565.2606  1    < 2e-16 ***
#> grade:sentence 349.7723  2    < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Because a latent model is no longer saturated, global fit measures become informative:

fit_measures(fit_latent)
#>          chisq             df         pvalue            cfi            tli 
#>   1.495041e+02   4.900000e+01   4.223177e-12   8.804705e-01   8.390011e-01 
#>          rmsea rmsea.ci.lower rmsea.ci.upper           srmr 
#>   1.104941e-01   9.053390e-02   1.309515e-01   7.270585e-02

Covariates

covariates adds covariates. Each one enters as a predictor of every contrast variable, so you get a main effect for the covariate plus its interactions with the within-subject terms. Those interactions say whether the experimental effects differ depending on the covariate.

fit_cov <- semnova(
    data = reading_two,
    id = "id",
    dv = "dv",
    indicator = "indicator",
    withins = c("grade", "sentence"),
    covariates = "iq",
    contrasts_arg = list(grade = "contr.sum", sentence = "contr.sum")
)

anova_table(fit_cov)
#> Hypothesis tests:
#>                     Chisq Df Pr(>Chisq)    
#> (Intercept)        1.4807  1    0.22367    
#> grade             29.0767  2  4.854e-07 ***
#> sentence          19.6252  1  9.422e-06 ***
#> grade:sentence    21.0598  2  2.673e-05 ***
#> iq                 2.3189  1    0.12781    
#> grade:iq           4.0542  2    0.13172    
#> sentence:iq        2.6564  1    0.10313    
#> grade:sentence:iq  6.0968  2    0.04744 *  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Covariates may be categorical. A factor is coded into contrast columns for you, using contrasts_arg or contr.treatment by default. Its columns stay one term, so a three-level covariate is a single effect with two degrees of freedom rather than two separate effects:

reading_grouped <- reading_manifest
reading_grouped$iq_group <- cut(reading_grouped$iq, 3, labels = c("low", "mid", "high"))

fit_catcov <- semnova(
    data = reading_grouped,
    id = "id",
    dv = "re_fix_time",
    withins = c("grade", "sentence"),
    covariates = "iq_group",
    contrasts_arg = list(grade = "contr.sum", sentence = "contr.sum")
)

anova_table(fit_catcov)
#> Hypothesis tests:
#>                           Chisq Df Pr(>Chisq)    
#> (Intercept)              2.6759  1  0.1018756    
#> grade                    6.3551  2  0.0416880 *  
#> sentence                13.2009  1  0.0002798 ***
#> grade:sentence          11.1570  2  0.0037783 ** 
#> iq_group                 2.2087  2  0.3314266    
#> grade:iq_group          22.0234  4  0.0001983 ***
#> sentence:iq_group        0.9646  2  0.6173500    
#> grade:sentence:iq_group  8.1860  4  0.0849986 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

expand_categorical_covariates() is the function that does the coding, should you want the columns themselves:

expand_categorical_covariates(reading_grouped, "iq_group", list())$terms
#> $iq_group
#> [1] "iq_group_mid"  "iq_group_high"

Measurement invariance

With a latent dependent variable you can ask whether the measurement model holds equally across measurement occasions. invariance_within takes the usual four levels, each nested in the previous one.

invariance_fit <- function(level) {
    semnova(
        data = reading_two,
        id = "id",
        dv = "dv",
        indicator = "indicator",
        withins = c("grade", "sentence"),
        invariance_within = level,
        contrasts_arg = list(grade = "contr.sum", sentence = "contr.sum")
    )
}

configural <- invariance_fit("configural")
weak <- invariance_fit("weak")
strong <- invariance_fit("strong")

Passing several fits to anova() compares them as nested models, which is how you test whether a stricter level of invariance is tenable. Alongside the chi-square difference test it reports each model’s own RMSEA and CFI and the change in each from the previous model. A difference test rejects almost any constraint in a large enough sample, so invariance is conventionally judged on those changes as well. indices takes any measure lavaan::fitMeasures() knows, in the order you want the columns: indices = c("rmsea", "cfi", "srmr"), or NULL for none. (lavaan’s RMSEA_D column is a third thing again: the RMSEA of the difference test itself.)

anova(configural, weak, strong)
#> 
#> Chi-Squared Difference Test
#> 
#>            Df    AIC    BIC   RMSEA RMSEA diff     CFI   CFI diff  Chisq
#> configural 39 2968.2 3127.6 0.10843            0.90839            116.03
#> weak       44 2970.8 3114.5 0.10695 -0.0014764 0.89944 -0.0089499 128.56
#> strong     49 2981.7 3109.8 0.11049  0.0035405 0.88047 -0.0189651 149.50
#>            Chisq diff Df diff  RMSEA_D Pr(>Chisq)    
#> configural                                           
#> weak           12.525       5 0.094651  0.0282572 *  
#> strong         20.947       5 0.137782  0.0008291 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

invariance_between does the same across groups in a multigroup model.

Identification is handled for you, including a case that is easy to get wrong: a construct measured at a single occasion by exactly two indicators needs both loadings fixed to 1 and both intercepts fixed to 0. Partial fixing (the usual referent-indicator or effects-coding scheme) leaves it under-identified. identification_scheme() is the function that decides:

identification_scheme(n_occasions = 1, n_indicators = 2, lv_scaling = "effect")
#> $method
#> [1] "full_tau_equivalence"
identification_scheme(n_occasions = 6, n_indicators = 2, lv_scaling = "effect")
#> $method
#> [1] "effect_coding"

Residual structure

With a latent dependent variable measured repeatedly, the same indicator’s residuals usually correlate over occasions, and their variances may be worth constraining. Both are described with resid_set(), which says which residual terms a constraint covers: the indicators, and the conditions.

resid_cov adds residual covariances among a set’s members. Name an indicator and it covers every occasion:

fit_rescov <- semnova(
    data = reading_two,
    id = "id",
    dv = "dv",
    indicator = "indicator",
    withins = c("grade", "sentence"),
    resid_cov = resid_set("fix_count"),
    contrasts_arg = list(grade = "contr.sum", sentence = "contr.sum")
)

anova(fit_latent, fit_rescov)
#> 
#> Chi-Squared Difference Test
#> 
#>            Df    AIC    BIC    RMSEA RMSEA diff     CFI  CFI diff   Chisq
#> fit_rescov 34 2913.4 3088.4 0.054931            0.97950            51.235
#> fit_latent 49 2981.7 3109.8 0.110494   0.055563 0.88047 -0.099031 149.504
#>            Chisq diff Df diff RMSEA_D Pr(>Chisq)    
#> fit_rescov                                          
#> fit_latent     98.269      15 0.18178  2.775e-14 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Covariances in a set are free by default; equal = TRUE makes them one parameter. Conditions restrict the set to part of the design. Here it covers residual covariances among landolt sentences only, constrained equal:

fit_landolt <- semnova(
    data = reading_two,
    id = "id",
    dv = "dv",
    indicator = "indicator",
    withins = c("grade", "sentence"),
    resid_cov = resid_set("fix_count", sentence = "landolt", equal = TRUE),
    contrasts_arg = list(grade = "contr.sum", sentence = "contr.sum")
)

fit_measures(fit_landolt)
#>          chisq             df         pvalue            cfi            tli 
#>   1.245263e+02   4.800000e+01   1.019855e-08   9.089873e-01   8.748575e-01 
#>          rmsea rmsea.ci.lower rmsea.ci.upper           srmr 
#>   9.741598e-02   7.659523e-02   1.185954e-01   6.414789e-02

resid_var uses the same sets to constrain residual variances to be equal. Omit the indicator and the set covers all of them; TRUE is shorthand for one set covering everything:

equal_over_time <- semnova(
    data = reading_two,
    id = "id",
    dv = "dv",
    indicator = "indicator",
    withins = c("grade", "sentence"),
    # each indicator's residual variance held constant across occasions
    resid_var = list(resid_set("fix_count"), resid_set("ini_land_pos")),
    contrasts_arg = list(grade = "contr.sum", sentence = "contr.sum")
)

anova(fit_latent, equal_over_time)
#> 
#> Chi-Squared Difference Test
#> 
#>                 Df    AIC    BIC   RMSEA RMSEA diff     CFI  CFI diff  Chisq
#> fit_latent      49 2981.7 3109.8 0.11049            0.88047           149.50
#> equal_over_time 59 3009.6 3106.5 0.11818  0.0076804 0.83537 -0.045097 197.42
#>                 Chisq diff Df diff RMSEA_D Pr(>Chisq)    
#> fit_latent                                               
#> equal_over_time     47.919      10 0.15024  6.421e-07 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A set is printed the way you wrote it, and is checked as you write it rather than deep inside the fit:

resid_set("fix_count", sentence = "landolt", equal = TRUE)
#> <resid_set> fix_count | sentence = landolt | equal
try(resid_set("fix_count", grade = "grade7"))
#> <resid_set> fix_count | grade = grade7 | equal or free, by argument

Two cases where a residual constraint cannot apply, both warned about rather than silently dropped: a manifest dependent variable has no residuals of its own (the observed variable is the construct), and a set covering fewer than two residual terms, which happens with a single occasion, has nothing to constrain.

build_resid_cov_syntax() and build_resid_var_syntax() write this part of the model, if you want the syntax itself.

cat(build_resid_var_syntax(list(c("Y1", "Y2", "Y3"))))
#> Y1 ~~ .theta_set_1_1*Y1
#> Y2 ~~ .theta_set_1_1*Y2
#> Y3 ~~ .theta_set_1_1*Y3

The covariance structure of the contrast variables

By default the contrast variables have an unstructured covariance matrix: every variance and covariance is freely estimated, and the model is saturated. Two options impose the classical repeated-measures structures instead.

sphericity constrains a model term’s contrast variables to equal variances and zero covariances among themselves. That is the assumption a classical RM-ANOVA makes, stated explicitly here and therefore testable. You name the terms, in the same vocabulary as the ANOVA table: TRUE for every term wide enough to constrain, or particular terms such as "grade" and "grade:sentence" (an interaction may also be given as its factors, c("grade", "sentence")).

fit_spher <- semnova(
    data = reading_manifest,
    id = "id",
    dv = "re_fix_time",
    withins = c("grade", "sentence"),
    sphericity = TRUE,
    contrasts_arg = list(grade = "contr.sum", sentence = "contr.sum")
)

fit_measures(fit_spher)
#>          chisq             df         pvalue            cfi            tli 
#>    125.6318847      4.0000000      0.0000000      0.6790877     -0.2034213 
#>          rmsea rmsea.ci.lower rmsea.ci.upper           srmr 
#>      0.4254407      0.3632950      0.4909841      0.2467095

Because the sphericity model is nested in the unstructured one, anova() tests the assumption directly:

anova(fit, fit_spher)
#> 
#> Chi-Squared Difference Test
#> 
#>           Df    AIC    BIC   RMSEA RMSEA diff     CFI CFI diff  Chisq
#> fit        0 1494.6 1579.0 0.00000            1.00000            0.00
#> fit_spher  4 1612.2 1684.1 0.42544    0.42544 0.67909 -0.32091 125.63
#>           Chisq diff Df diff RMSEA_D Pr(>Chisq)    
#> fit                                                
#> fit_spher     125.63       4 0.42544  < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

compound_symmetry puts the structure on the occasions instead: one variance for every occasion and one covariance for every pair, with the contrast variables’ own (co)variances fixed to zero.

fit_cs <- semnova(
    data = reading_manifest,
    id = "id",
    dv = "re_fix_time",
    withins = c("grade", "sentence"),
    compound_symmetry = TRUE,
    contrasts_arg = list(grade = "contr.sum", sentence = "contr.sum")
)

anova(fit, fit_cs)
#> 
#> Chi-Squared Difference Test
#> 
#>        Df    AIC    BIC   RMSEA RMSEA diff     CFI CFI diff  Chisq Chisq diff
#> fit     0 1494.6 1579.0 0.00000            1.00000            0.00           
#> fit_cs 19 1730.6 1755.6 0.28263    0.28263 0.32728 -0.67272 273.97     273.97
#>        Df diff RMSEA_D Pr(>Chisq)    
#> fit                                  
#> fit_cs      19 0.28263  < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

In a multigroup model, variance_homogeneity = TRUE constrains those (co)variances to be equal across groups, and group_weights = "stochastic" treats group membership as random rather than fixed by design. Each group’s frequency is then estimated and reported as a relfreq definition.

fit_weighted <- semnova(
    data = reading_manifest,
    id = "id",
    dv = "re_fix_time",
    withins = c("grade", "sentence"),
    betweens = "sex",
    group_weights = "stochastic",
    variance_homogeneity = TRUE,
    contrasts_arg = list(
        grade = "contr.sum", sentence = "contr.sum", sex = "contr.treatment"
    )
)

estimates <- lavaan::parTable(fit_weighted$sem_obj)
estimates[estimates$op == ":=", c("lhs", "est")]
#>          lhs     est
#> 183        N 168.000
#> 184 relfreq1   0.536
#> 185 relfreq2   0.464

build_pi_layer_syntax() and build_group_weight_syntax() are the functions that write this part of the model, should you want the syntax itself.

cat(build_group_weight_syntax(c("m", "w"), "stochastic"))
#> group % c(.gw_1, .gw_2)*w
#> N := exp(.gw_1) + exp(.gw_2)
#> relfreq1 := exp(.gw_1) / N
#> relfreq2 := exp(.gw_2) / N

Hypothesis tests

Every row of the ANOVA table is a hypothesis, and every hypothesis can be inspected on its own. test_hypothesis() returns both the individual contrast estimates and the omnibus test.

term_names <- vapply(fit$hypotheses, function(h) h$description, character(1))
term_names
#> [1] "(Intercept)"    "grade"          "sentence"       "grade:sentence"

grade_effect <- fit$hypotheses[[which(term_names == "grade")]]
test_hypothesis(fit, grade_effect)
#> Hypothesis: grade
#> 
#> Contrast estimates:
#>        estimate std_error z value ci_lower ci_upper sd_between Pr(>|z|)    
#> grade1   0.5915   0.02330  25.384   0.5458   0.6372      0.302  < 2e-16 ***
#> grade2  -0.0557   0.01636  -3.404  -0.0877  -0.0236      0.212 0.000664 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Omnibus Wald test: chisq(2) = 832.066, p = 0.0000

contrast_estimates() is the same function under a name that describes the other half of what it returns; use whichever reads better at the call site. sd_between is the between-person standard deviation of the effect. It says how much the effect varies across people, not how precisely it is estimated.

contrast_estimates(fit, grade_effect, level = 0.99)$table
#> Contrast estimates:
#>        estimate std_error z value ci_lower ci_upper sd_between  Pr(>|z|)    
#> grade1  0.59151  0.023302  25.384  0.53148  0.65153    0.30203 < 2.2e-16 ***
#> grade2 -0.05567  0.016355  -3.404 -0.09780 -0.01355    0.21199  0.000664 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

contrast_table() stacks the estimates for every term at once.

contrast_table(fit)
#> Contrast estimates:
#>                                  estimate std_error z value ci_lower ci_upper
#> (Intercept)                      -0.06773  0.031207 -2.1704 -0.12890 -0.00657
#> grade: grade1                     0.59151  0.023302 25.3840  0.54583  0.63718
#> grade: grade2                    -0.05567  0.016355 -3.4040 -0.08773 -0.02362
#> sentence: sentence1               0.54444  0.023983 22.7006  0.49743  0.59144
#> grade:sentence: grade1:sentence1  0.27811  0.019907 13.9709  0.23910  0.31713
#> grade:sentence: grade2:sentence1 -0.02717  0.014407 -1.8861 -0.05541  0.00106
#>                                  sd_between  Pr(>|z|)    
#> (Intercept)                         0.40449  0.029976 *  
#> grade: grade1                       0.30203 < 2.2e-16 ***
#> grade: grade2                       0.21198  0.000664 ***
#> sentence: sentence1                 0.31086 < 2.2e-16 ***
#> grade:sentence: grade1:sentence1    0.25802 < 2.2e-16 ***
#> grade:sentence: grade2:sentence1    0.18673  0.059276 .  
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Custom hypotheses

You are not limited to the terms of the design. Build any general linear hypothesis with hypothesis() and pass it to the same function. There is no separate code path for custom tests.

A hypothesis is specified by up to four matrices: M picks out within-subject contrasts, G group contrasts, L covariate contrasts, and K the value tested against (zero by default). Any of them may be omitted when the design makes them unnecessary.

Here we test a single within-subject contrast on its own:

M <- matrix(0, nrow = length(fit$pis), ncol = 1)
M[2, 1] <- 1

my_hypothesis <- hypothesis(M = M, description = "second contrast only")
test_hypothesis(fit, my_hypothesis)
#> Hypothesis: second contrast only
#> 
#> Contrast estimates:
#>           estimate std_error z value ci_lower ci_upper sd_between Pr(>|z|)    
#> contrast1   0.5915    0.0233   25.38   0.5458   0.6372      0.302   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Omnibus Wald test: chisq(1) = 644.345, p = 0.0000

validate_hypothesis() checks a hypothesis against a fit and fills in whatever can be defaulted, which is useful when building hypotheses programmatically:

validate_hypothesis(my_hypothesis, fit)
dim(my_hypothesis$K)
#> [1] 1 1 1

Plots

plot() visualises one term’s contrast estimates. Name the term with terms= (order does not matter) or pass a hypothesis object directly.

plot(fit, terms = c("grade", "sentence"), interval = "ci")

Switching interval to "sd" shows between-person variability instead of the confidence interval. It shows the spread of the effect across people.

plot(fit, terms = "grade", interval = "sd", geom = c("point", "line"))

dag() draws the model as a SEM path diagram, laid out in the structure a semnova model has: residual terms, then indicators, then the latent variables they measure, then the contrast variables, then any covariates.

Standard SEM conventions apply. Single-headed arrows are directed effects (loadings, regressions); double-headed arrows and self-loops are (co)variances. Squares are manifest variables, circles latent ones, triangles the constant. Fixed parameters are marked with a ring accent over the estimate.

By default the diagram shows the measurement and regression structure. That is the part you usually want to see. A fully crossed design has a great many (co)variances, so drawing everything at once is rarely readable.

dag(fit_cov)

what controls which parts are drawn. Give any of "measurement", "regression", "variance", "covariance", "resid_var", "resid_cov", "intercept", or "all". Residual variances appear as their own epsilon nodes, and intercepts as triangles:

dag(fit_cov, what = c("measurement", "regression", "resid_var", "variance"))

Multigroup fits are drawn as facets, one panel per group, labelled with the group name:

dag(fit_between)

labels controls the estimates: "all" (the default) labels every parameter drawn. Fixed ones, such as the contrast loadings, are shown inside a ring accent, so design and estimate stay distinguishable. "free" labels only freely estimated parameters, which thins out a dense diagram; "none" gives a bare structural one.

If a dense model looks cluttered, x_scale and y_scale stretch the node positions horizontally and vertically. Only their ratio matters, because stretching both equally gives the same diagram. The tighter direction keeps its spacing and the other grows away from it. Here the columns are pushed apart while the rows stay as they were:

p <- dag(fit_cov, x_scale = 1.5)
attr(p, "fig_size") # the figure size this needs, in inches
#>  width height 
#>   12.9    8.1
p

A scaled diagram is drawn at a physical size rather than fitted to the figure, so it needs a figure big enough to hold it. That is what the fig_size attribute is for, and the chunk options above were set from it. Unscaled, dag() fits whatever figure it is given, and a bigger fig.width/fig.height is the simplest cure for crowding.

The curved paths, meaning variance self-loops and covariance arcs, are drawn at a fixed size and a fixed bend angle. so scaling the layout or resizing the figure leaves them alone. curve_scale is what changes them, with variance_scale and covariance_scale following it unless you set them yourself:

dag(fit_between, what = c("measurement", "regression", "variance", "covariance"),
    curve_scale = 1.6)

A scaled diagram is fitted to the figure rather than cropped by it: if the figure is smaller than fig_size, you get a smaller diagram, not a piece of one. Since the markers and text keep their fixed size, though, a figure well under fig_size will crowd again. Scaling sets the shape, the figure gives the room.

dag_save() writes a diagram to a file and crops the whitespace off it. There is always some: the layout is fitted to the figure, and the margins have to be wide enough for the labels that overhang the panel. Dimensions default to the plot’s own fig_size.

cropped <- tempfile(fileext = ".png")
plain <- tempfile(fileext = ".png")

dag_save(cropped, dag(fit_cov), dpi = 100)
dag_save(plain, dag(fit_cov), dpi = 100, trim = FALSE)

vapply(c(cropped = cropped, uncropped = plain), file.size, numeric(1))
#>   cropped uncropped 
#>    138006    140502

Cropping needs the magick package and a raster format; without either the file is still written, uncropped, with a message saying why. padding sets the margin left around the diagram, in pixels.

Looking inside a fit

model_parameters() is the single canonical read of a fitted model: every parameter, classified by type.

params <- model_parameters(fit_cov)
table(params$type)
#> 
#>  constraint  covariance   intercept measurement  regression   resid_var 
#>          22          30          25          48           6          13 
#>    variance 
#>          12

head(model_parameters(fit_cov, what = "regression")[c("lhs", "rhs", "est", "se", "type")])
#>      lhs rhs           est           se       type
#> 127 .pi1  iq -0.0036880775 0.0024219231 regression
#> 128 .pi2  iq -0.0030748545 0.0016105930 regression
#> 129 .pi3  iq  0.0004591676 0.0010789739 regression
#> 130 .pi4  iq -0.0026161540 0.0016051477 regression
#> 131 .pi5  iq -0.0022663795 0.0012808644 regression
#> 132 .pi6  iq -0.0004644074 0.0009207377 regression

regression_coefficients() reshapes those into the arrays the contrast machinery uses: the coefficients, their sampling covariance, and the between-person covariance of the contrast variables.

coefs <- regression_coefficients(fit_cov)
dim(coefs$B_array)
#> [1] 2 6 1
dimnames(coefs$B_array)[["covariate"]]
#> [1] "(Intercept)" "iq"

The underlying lavaan fit is always available if you need something the package does not expose:

class(fit_cov$sem_obj)
#> [1] "lavaan"
#> attr(,"package")
#> [1] "lavaan"

The lower-level interface

semnova() derives the latent variables and contrast matrix from your design. When you want to specify them yourself, for a contrast matrix no standard contrast function produces, use lgc() directly.

You supply data in wide format, the constructs and their indicators, and the contrast matrix.

wide <- fit$data
custom_contrasts <- solve(cbind(1, stats::contr.poly(6)))

fit_lgc <- lgc(
    data = wide,
    etas = fit$etas_mmodel,
    c_matrix_within = custom_contrasts,
    hypotheses = hypothesis(M = M, description = "custom contrast")
)

anova_table(fit_lgc)
#> Hypothesis tests:
#>                  Chisq Df Pr(>Chisq)    
#> custom contrast 876.75  1  < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Everything shown above works the same way on an lgc() fit. The ANOVA-style interface is a convenience layer, not a separate model.

Building blocks

The syntax generators semnova() calls are exported too, for anyone who wants to inspect or reuse the generated lavaan syntax. build_measurement_syntax() writes the measurement model:

etas <- list(eta1 = c("Y1", "Y2", "Y3"))
cat(build_measurement_syntax(etas, lv_scaling = "effect"))
#> ## LOADINGS:
#> 
#> eta1 =~ NA*Y1 + .lambda_1_1_1*Y1
#> eta1 =~ NA*Y2 + .lambda_1_2_1*Y2
#> eta1 =~ NA*Y3 + .lambda_1_3_1*Y3
#> 
#> ## INTERCEPTS:
#> 
#> Y1 ~ NA*1 + .nu_1_1_1*1
#> Y2 ~ NA*1 + .nu_1_2_1*1
#> Y3 ~ NA*1 + .nu_1_3_1*1
#> 
#> ## UNIQUE-FACTOR VARIANCES:
#> 
#> Y1 ~~ NA*Y1 + .theta_1_1_1*Y1
#> Y2 ~~ NA*Y2 + .theta_1_2_1*Y2
#> Y3 ~~ NA*Y3 + .theta_1_3_1*Y3
#> 
#> ## LATENT MEANS/INTERCEPTS:
#> 
#> eta1 ~ NA*1 + .alpha_1_1*1
#> 
#> ## COMMON-FACTOR VARIANCES:
#> 
#> eta1 ~~ NA*eta1 + .psi_1_1*eta1
#> 
#> ## MODEL CONSTRAINTS:
#> 
#> .lambda_1_1_1 == 3 - .lambda_1_2_1 - .lambda_1_3_1
#> .nu_1_1_1 == 0 - .nu_1_2_1 - .nu_1_3_1

build_pi_layer_syntax() writes the contrast layer, build_covariate_syntax() the covariate model and its regressions, and build_resid_cov_syntax() any residual covariances you want to allow:

cat(build_pi_layer_syntax(matrix(c(1, 1, 1, -1), nrow = 2), eta_names = c("eta1", "eta2")))
#> .pi1 =~ 0.5*eta1 + 0.5*eta2
#> .pi2 =~ 0.5*eta1 + -0.5*eta2
#> .pi1 ~ .alpha_pi_1_1*1
#> .pi2 ~ .alpha_pi_2_1*1
#> .pi1 ~~ .psi_pi_1_1*.pi1
#> .pi1 ~~ .psi_pi_1_2_1*.pi2
#> .pi2 ~~ .psi_pi_2_1*.pi2
#> eta1 ~ 0*1
#> eta2 ~ 0*1
#> eta1 ~~ 0*eta1
#> eta1 ~~ 0*eta2
#> eta2 ~~ 0*eta2

cat(build_resid_cov_syntax(equal_resid_cov = list(c("Y1", "Y2", "Y3"))))
#> Y1 ~~ .theta_cov_1_1*Y2
#> Y1 ~~ .theta_cov_1_1*Y3
#> Y2 ~~ .theta_cov_1_1*Y3

strip_construct_structure() removes a construct’s own mean and variance lines, which is what lets the measurement model and the contrast layer be combined without specifying the same parameter twice:

cat(strip_construct_structure("eta1 =~ Y1\neta1 ~~ eta1\neta1 ~ 1", "eta1"))
#> eta1 =~ Y1

parse_design() turns factor columns into the contrast matrices a design implies, and generate_anova_hypotheses() crosses two designs into one hypothesis per model term:

within_design <- parse_design(
    reading_manifest,
    factors = c("grade", "sentence"),
    contrasts_arg = list(grade = "contr.sum", sentence = "contr.sum")
)
within_design$effect_names
#> [1] "(Intercept)"    "grade"          "sentence"       "grade:sentence"
round(within_design$c_matrix, 3)
#>                       1      2      3      4      5      6
#> (Intercept)       0.167  0.167  0.167  0.167  0.167  0.167
#> grade1            0.333 -0.167 -0.167  0.333 -0.167 -0.167
#> grade2           -0.167  0.333 -0.167 -0.167  0.333 -0.167
#> sentence1         0.167  0.167  0.167 -0.167 -0.167 -0.167
#> grade1:sentence1  0.333 -0.167 -0.167 -0.333  0.167  0.167
#> grade2:sentence1 -0.167  0.333 -0.167  0.167 -0.333  0.167

contrast_matrix() builds the Kronecker contrast matrix behind every test:

R <- contrast_matrix(G = matrix(1), L = matrix(1), M = M)
dim(R$R)
#> [1] 1 6

Parameter names

Parameters carry structured labels, generated and read back by one small set of functions so that writing the model and reading its results can never disagree.

A parameter’s identity is its type and its indices, and every rendering comes from those: the lavaan label the model is written with, a readable text form, and the plotmath dag() draws on the arrows.

param_label("regression", pi = 2, covariate = 1, group = 1)
#> [1] ".beta_2_1_1"
param_label("regression", pi = 2, covariate = 1, group = 1, format = "text")
#> [1] "beta_{2,1,1}"
param_label("regression", pi = 2, covariate = 1, group = 1, format = "plotmath")
#> [1] "beta[2*','*1*','*1]"
parse_param_label(".beta_2_1_1")
#> $type
#> [1] "regression"
#> 
#> $pi
#> [1] 2
#> 
#> $covariate
#> [1] 1
#> 
#> $group
#> [1] 1

The labels also state their equality constraints. Every parameter of the measurement model is labelled by the same engine. An index of 0, shown as *, means “one parameter across this whole dimension”. Under strong within-invariance the six occasions share one loading per indicator, and the label says so:

labels <- lavaan::parTable(fit_latent$sem_obj)$label
loadings <- unique(grep("^\\.lambda_", labels, value = TRUE))
loadings
#> [1] ".lambda_0_1_1" ".lambda_0_2_1"
pretty_label(loadings[1], format = "text") # * = shared across occasions
#> [1] "lambda_{*,1,1}"

pretty_label() is for labels you already have as a string, such as those read off a fitted model; prefer param_label(format = ) whenever you know the parameter’s identity.

regression_labels() enumerates one family as a positioned array: the pi-layer regression coefficients, (covariate x pi x group). That is the array regression_coefficients() fills with estimates, and the two come from the same call so they cannot drift apart. Everything else in the model is named one label at a time, where it is written.

regression_labels(pis = 1:2, covariates = "iq", groups = "1")[, , 1]
#>              pi
#> covariate     1               2              
#>   (Intercept) ".alpha_pi_1_1" ".alpha_pi_2_1"
#>   iq          ".beta_1_1_1"   ".beta_2_1_1"

n_occasions() reports how many occasions a construct is measured at. The identification rules depend on that count.

n_occasions("iq", n_within_occasions = 6, covariate_names = "iq")
#> [1] 1

What is not here yet

These are known gaps rather than oversights:

  • Power analysis. Deferred. It will return in a later release.
  • Marginal effects. The old Effect() method is not carried over. Use plot() on a hypothesis instead.

Session information

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] semnova_0.3.0.9001 rmarkdown_2.32    
#> 
#> loaded via a namespace (and not attached):
#>  [1] viridis_0.6.5      sass_0.4.10        generics_0.1.4     tidyr_1.3.2       
#>  [5] digest_0.6.39      magrittr_2.0.5     evaluate_1.0.5     grid_4.6.1        
#>  [9] RColorBrewer_1.1-3 fastmap_1.2.0      jsonlite_2.0.0     ggrepel_0.9.8     
#> [13] gridExtra_2.3.1    purrr_1.2.2        viridisLite_0.4.3  scales_1.4.0      
#> [17] tweenr_2.0.3       pbivnorm_0.6.0     jquerylib_0.1.4    mnormt_2.1.2      
#> [21] cli_3.6.6          graphlayouts_1.2.5 rlang_1.3.0        polyclip_1.10-7   
#> [25] tidygraph_1.3.1    withr_3.0.3        cachem_1.1.0       yaml_2.3.12       
#> [29] otel_0.2.0         tools_4.6.1        parallel_4.6.1     memoise_2.0.1     
#> [33] dplyr_1.2.1        ggplot2_4.0.3      buildtools_1.0.0   vctrs_0.7.3       
#> [37] R6_2.6.1           magick_2.9.1       stats4_4.6.1       lifecycle_1.0.5   
#> [41] MASS_7.3-66        ggraph_2.2.2       pkgconfig_2.0.3    bslib_0.12.0      
#> [45] pillar_1.11.1      gtable_0.3.6       Rcpp_1.1.2         glue_1.8.1        
#> [49] ggforce_0.5.0      xfun_0.62          tibble_3.3.1       tidyselect_1.2.1  
#> [53] sys_3.4.3          knitr_1.52         farver_2.1.2       htmltools_0.5.9   
#> [57] igraph_2.3.4       lavaan_0.7-3       maketools_1.3.2    labeling_0.4.3    
#> [61] compiler_4.6.1     S7_0.2.2           quadprog_1.5-8     semTools_0.5-10