Nonparametric Longitudinal Data Analysis with nparLD

library(nparLD)

Introduction

The nparLD package provides nonparametric methods for longitudinal and repeated-measures data in factorial experiments. It implements inference for hypotheses in marginal distribution functions and in unweighted relative marginal effects. The redesigned interface uses a general formula syntax and supports crossed factorial designs, incomplete observations, dependent replicates, multiple contrast procedures, and simultaneous confidence intervals.

The main function is nparLD(). It supports crossed factorial designs with whole-plot and subplot factors, incomplete observations, dependent replicate measurements, rank- and pseudo-rank-based inference, tests in marginal distribution functions, tests in unweighted relative marginal effects, and graphical summaries of estimated effects and contrasts.

Data structure and formula interface

The data must be supplied in long format. Each row corresponds to one observed measurement. At minimum, the data frame contains a response variable, a subject identifier, and one or more design factors. Repeated measurements are represented by one or more within-subject factors. Whole-plot factors, such as treatment groups or strata, vary between subjects.

A typical call has the form


nparLD(response ~ factor1 * factor2,
       data = dat,
       subject = "subject")

The formula determines the factorial structure of the analysis. Interactions are specified in the usual R formula notation. Factor ordering is controlled by the levels of the corresponding variables in the input data. Thus, when a specific ordering of time points, treatment groups, or ordinal categories is required, the variables should be converted to factors with the desired level order before calling nparLD().

Classical designs and the formula interface

Earlier versions of nparLD used design-specific functions for common longitudinal layouts. In the redesigned version, these designs are specified through the same formula interface allowing for missing values and possibly dependent replicates.

Classical design Formula
LD-F1 resp ~ time
F1-LD-F1 resp ~ group * time
LD-F2 resp ~ time1 * time2
F1-LD-F2 resp ~ group * time1 * time2
F2-LD-F1 resp ~ group1 * group2 * time

Factor ordering is controlled by setting factor levels in the input data before calling nparLD().

Marginal distributions and relative effects

For each factorial cell, the basic distributional object is the marginal distribution of the response in that cell. The package uses the normalized marginal distribution function

\[ \frac{1}{2}{F_{is}^{+}(x) + F_{is}^{-}(x)}, \]

where

\[ F_{is}^{+}(x) = P(X_{is} \le x), \qquad F_{is}^{-}(x) = P(X_{is} < x). \]

This convention corresponds to mid-ranks and provides a common treatment of continuous, discrete, ordinal, and tied responses.

Relative marginal Mann-Whitney effects describe the position of a marginal distribution relative to a reference distribution. Values larger than 0.5 indicate a tendency toward larger responses than the reference, while values smaller than 0.5 indicate a tendency toward smaller responses. In factorial longitudinal designs, the choice of the reference distribution is important and leads to the distinction between weighted and unweighted effects.

Weighted and unweighted relative effects

Relative effects are defined with respect to a reference distribution. In a factorial longitudinal design, this reference distribution can be formed in two different ways.

The weighted relative effect uses the sample-size weighted reference distribution. Marginal distributions from cells with larger sample sizes therefore contribute more strongly to the reference distribution than cells with smaller sample sizes. Weighted effects are closely connected to ordinary ranks: they describe the relative position of a cell distribution within the empirical pooled distribution induced by the observed sample sizes.

This makes weighted effects useful as descriptive summaries of the observed data. However, they depend on the allocation of sample sizes across cells. In an unbalanced design, the same set of marginal distributions could lead to different weighted effects if the sample-size allocation were changed. For this reason, weighted effects are not used as the target of hypotheses in relative marginal effects.

The unweighted relative effect uses an equally weighted reference distribution. Each factorial cell contributes equally to the reference distribution, independently of its sample size. The resulting effects therefore refer to the factorial marginal distributions themselves rather than to the particular sample-size allocation of the study. This is the effect scale used for hypotheses in relative marginal effects.

In nparLD, hypotheses specified by

hypothesis = "H0p"

are hypotheses in unweighted relative marginal effects. These effects are estimated by pseudo-ranks. Weighted effects, in contrast, are reported as descriptive rank-based summaries and are naturally associated with hypotheses in marginal distribution functions,

hypothesis = "H0F"

when the interest lies in comparing the underlying marginal distributions.

The distinction is most important in unbalanced designs. In balanced designs, weighted and unweighted reference distributions often lead to similar interpretations. In unbalanced designs, however, weighted effects partly reflect the empirical sample-size allocation, whereas unweighted effects keep the factorial cells on equal footing.

Hypotheses H0F and H0p

The argument hypothesis = “H0F” specifies hypotheses in marginal distribution functions. In abstract notation, these hypotheses have the form

\[ H_0^F : \mathbf{C F} = \mathbf{0}, \]

where \(\mathbf{F}\) collects the marginal distribution functions and \(\mathbf{C}\) is a contrast matrix. Such hypotheses compare complete marginal distributions.

The argument hypothesis = “H0p” specifies hypotheses in unweighted relative marginal effects,

\[ H_0^p : \mathbf{C \psi} = \mathbf{0}, \]

where \(\mathbf{\psi}\) collects the unweighted relative marginal effects. These hypotheses focus on the relative positions of the marginal distributions with respect to a common unweighted reference distribution. They are especially useful for effect interpretation, multiple contrasts, simultaneous confidence intervals, and plots. Note that this null hypothesis does not imply equal marginal distributions and allows for unequal variances and/or other moments even under the null hypothesis. Therefore, testing "H0p" is also known as the nonparametric Behrens-Fisher problem.

The two hypothesis types therefore answer related but distinct questions. H0F asks whether marginal distribution functions are equal in the sense specified by the contrast matrix. H0p asks whether the corresponding unweighted relative marginal effects satisfy the specified contrast.

A one-factor longitudinal example

data(dental)

fit_dental <- nparLD(
  resp ~ time,
  data = dental,
  subject = "subject",
  hypothesis = "H0p",
  covariance = TRUE
)

fit_dental
#> 
#> Nonparametric Longitudinal Analysis
#> -------------------------------------------------
#> Design: LD-F1-Design 
#> 
#> Subplot factors: time
#> Number of subjects: 16
#> Hypotheses: Tested in marginal Mann-Whitney effects (relative marginal effects)
#> Ranking: Pseudo-ranks
#> Confidence interval method: logit
#> 
#> Estimated relative effects
#> -------------------------------------------------
#>   time Nsubj lambda Nobs Nmiss ScoreMean    RTE std.error  lower  upper
#> 2   10    16     16   16     0     24.31 0.3721    0.1331 0.3095 0.4392
#> 3   12    16     16   16     0     37.03 0.5708    0.1379 0.5023 0.6367
#> 4   14    16     16   16     0     49.59 0.7671    0.1239 0.7010 0.8223
#> 1    8    16     16   16     0     19.06 0.2900    0.1330 0.2294 0.3592
#> 
#> Global test results
#> -------------------------------------------------
#> 
#> Wald-type statistic (WTS)
#>      Statistic df p-value
#> time  203.6774  3       0
#> 
#> ANOVA-type statistic (ATS)
#>      Statistic     df p-value
#> time    31.381 2.5496       0
#> 
#> Estimated covariance matrix
#> -------------------------------------------------
#>          1:8    1:10    1:12    1:14
#> 1:8   0.0177 -0.0053  0.0001 -0.0125
#> 1:10 -0.0053  0.0190 -0.0120 -0.0017
#> 1:12  0.0001 -0.0120  0.0154 -0.0035
#> 1:14 -0.0125 -0.0017 -0.0035  0.0177

A factorial longitudinal example

data(shoulder)

fit_shoulder <- nparLD(
  resp ~ group1 * group2 * time,
  data = shoulder,
  subject = "subject",
  hypothesis = "H0p",
  contrast = list("group1:time")
)

fit_shoulder
#> 
#> Nonparametric Longitudinal Analysis
#> -------------------------------------------------
#> Design: F2-LD-F1-Design 
#> 
#> Whole-plot factors: group1, group2
#> Subplot factors: time
#> Number of subjects: 41
#> Hypotheses: Tested in marginal Mann-Whitney effects (relative marginal effects)
#> Ranking: Pseudo-ranks
#> Confidence interval method: logit
#> 
#> Estimated relative effects
#> -------------------------------------------------
#>    group1 group2 time Nsubj lambda Nobs Nmiss ScoreMean    RTE std.error  lower
#> 1       N      F    1    11     11   11     0    152.71 0.6187    0.5560 0.4410
#> 2       N      F    2    11     11   11     0    172.57 0.6995    0.4518 0.5466
#> 3       N      F    3    11     11   11     0    160.56 0.6507    0.4547 0.5024
#> 4       N      F    4    11     11   11     0    180.86 0.7332    0.4278 0.5845
#> 5       N      F    5    11     11   11     0    144.72 0.5863    0.5004 0.4297
#> 6       N      F    6    11     11   11     0    130.89 0.5301    0.4050 0.4068
#> 7       N      M    1     8      8    8     0    124.69 0.5048    0.5573 0.3400
#> 8       N      M    2     8      8    8     0    166.79 0.6760    0.5067 0.5068
#> 9       N      M    3     8      8    8     0    175.24 0.7103    0.5629 0.5149
#> 10      N      M    4     8      8    8     0    171.01 0.6931    0.5320 0.5123
#> 11      N      M    5     8      8    8     0    148.38 0.6011    0.5736 0.4202
#> 12      N      M    6     8      8    8     0    120.46 0.4877    0.4606 0.3512
#> 13      Y      F    1    14     14   14     0    121.82 0.4932    0.4350 0.3636
#> 14      Y      F    2    14     14   14     0     98.47 0.3983    0.3976 0.2848
#> 15      Y      F    3    14     14   14     0     87.35 0.3530    0.3388 0.2573
#> 16      Y      F    4    14     14   14     0     99.22 0.4013    0.3633 0.2967
#> 17      Y      F    5    14     14   14     0     70.29 0.2837    0.2157 0.2225
#> 18      Y      F    6    14     14   14     0     82.16 0.3319    0.2970 0.2480
#> 19      Y      M    1     8      8    8     0    105.22 0.4257    0.5490 0.2715
#> 20      Y      M    2     8      8    8     0    110.75 0.4482    0.5027 0.3036
#> 21      Y      M    3     8      8    8     0     85.13 0.3440    0.3695 0.2411
#> 22      Y      M    4     8      8    8     0     74.74 0.3018    0.2831 0.2227
#> 23      Y      M    5     8      8    8     0     89.98 0.3637    0.4337 0.2437
#> 24      Y      M    6     8      8    8     0     89.98 0.3637    0.4337 0.2437
#>     upper
#> 1  0.7695
#> 2  0.8180
#> 3  0.7746
#> 4  0.8429
#> 5  0.7271
#> 6  0.6498
#> 7  0.6686
#> 8  0.8090
#> 9  0.8500
#> 10 0.8293
#> 11 0.7581
#> 12 0.6260
#> 13 0.6237
#> 14 0.5238
#> 15 0.4622
#> 16 0.5157
#> 17 0.3541
#> 18 0.4281
#> 19 0.5958
#> 20 0.6021
#> 21 0.4640
#> 22 0.3947
#> 23 0.5036
#> 24 0.5036
#> 
#> Global test results
#> -------------------------------------------------
#> 
#> Wald-type statistic (WTS)
#>                    Statistic df p-value
#> group1               16.9539  1  0.0000
#> group2                0.0476  1  0.8274
#> time                 14.7195  5  0.0116
#> group1:group2         0.0318  1  0.8584
#> group1:time          22.6602  5  0.0004
#> group2:time          12.1939  5  0.0322
#> group1:group2:time    5.2195  5  0.3897
#> 
#> ANOVA-type statistic (ATS)
#>                    Statistic     df p-value
#> group1               16.9539 1.0000  0.0000
#> group2                0.0476 1.0000  0.8274
#> time                  3.3876 2.7526  0.0203
#> group1:group2         0.0318 1.0000  0.8584
#> group1:time           3.7277 2.7325  0.0133
#> group2:time           1.1548 2.7092  0.3234
#> group1:group2:time    0.4311 2.6790  0.7084
#> 
#> Multiple contrast test procedure
#> -------------------------------------------------
#> Factor or interaction: group1:time
#> SCI method: fisher
#> 
#> Global result
#> -------------------------------------------------
#>   Statistic p.value
#> 1    3.5025  0.0102
#> 
#> Local contrast results
#> -------------------------------------------------
#>      Estimate Std.Error   lower   upper Statistic p.value df
#> C 1   -0.0731    0.2123 -0.1641  0.0191   -2.2089  0.1550 21
#> C 2    0.0080    0.1538 -0.0589  0.0748    0.3317  0.9970 21
#> C 3    0.0417    0.1233 -0.0119  0.0951    2.1667  0.1677 21
#> C 4    0.0565    0.1034  0.0116  0.1012    3.5025  0.0102 21
#> C 5    0.0107    0.1436 -0.0517  0.0730    0.4772  0.9850 21
#> C 6   -0.0438    0.1467 -0.1072  0.0200   -1.9124  0.2593 21
#> C 7    0.0731    0.2123 -0.0191  0.1641    2.2089  0.1553 21
#> C 8   -0.0080    0.1538 -0.0748  0.0589   -0.3317  0.9970 21
#> C 9   -0.0417    0.1233 -0.0951  0.0119   -2.1667  0.1676 21
#> C 10  -0.0565    0.1034 -0.1012 -0.0116   -3.5025  0.0106 21
#> C 11  -0.0107    0.1436 -0.0730  0.0517   -0.4772  0.9850 21
#> C 12   0.0438    0.1467 -0.0200  0.1072    1.9124  0.2592 21
#> 
#> Contrast matrix: available in component `Contrast.Matrix`.
#> Use `print(x, show.matrix = TRUE)` to display it.
plot(fit_shoulder)

plot(fit_shoulder$MCTP)

Missing observations

Longitudinal studies often contain incomplete observation patterns. A subject may miss one or more visits, or a measurement may be unavailable for a particular condition. The redesigned nparLD() function allows such incomplete subject-condition cells. Missing response values should be coded in the data set in the usual R way, for example as NA.

If a subject has a missing response at a particular time point or condition, the corresponding row should still contain the relevant factor-level information. Thus, the subject identifier, time point, group, and any other design factors should remain available in the data set, while only the response value is set to NA. This allows nparLD() to recognize the intended longitudinal design and the incomplete subject-condition cell automatically.

Observed measurements contribute to the estimation of the corresponding marginal distributions and effects, while missing subject-condition cells simply do not contribute observations for that cell. The analysis is still based on the independent sampling units, that is, on subjects. Thus, incomplete subjects are not automatically removed from the data as long as they contribute observations to at least part of the design. This is useful in longitudinal applications where complete-case analysis would discard a substantial amount of information. The implemented inference procedures are valid under the missing completely at random mechanism.

set.seed(123)

dat_miss <- dental
dat_miss$resp[c(2, 7, 12)] <- NA

fit_miss <- nparLD(
  resp ~ time,
  data = dat_miss,
  subject = "subject",
  hypothesis = "H0p"
)

fit_miss
#> 
#> Nonparametric Longitudinal Analysis
#> -------------------------------------------------
#> Design: LD-F1-Design 
#> 
#> Subplot factors: time
#> Number of subjects: 16
#> Hypotheses: Tested in marginal Mann-Whitney effects (relative marginal effects)
#> Ranking: Pseudo-ranks
#> Confidence interval method: logit
#> 
#> Estimated relative effects
#> -------------------------------------------------
#>   time Nsubj lambda Nobs Nmiss ScoreMean    RTE std.error  lower  upper
#> 2   10    16     16   16     0     22.34 0.3580    0.1583 0.2846 0.4387
#> 3   12    16     16   16     0     34.63 0.5595    0.1444 0.4881 0.6286
#> 4   14    16     16   16     0     47.07 0.7634    0.1387 0.6889 0.8246
#> 1    8    16     13   13     3     19.97 0.3191    0.1363 0.2563 0.3893
#> 
#> Global test results
#> -------------------------------------------------
#> 
#> Wald-type statistic (WTS)
#>      Statistic df p-value
#> time   98.9615  3       0
#> 
#> ANOVA-type statistic (ATS)
#>      Statistic   df p-value
#> time    24.039 2.84       0

The descriptive table in the output lists the number of subjects and missing values for every condition.

Dependent replicate measurements

In some longitudinal experiments, more than one measurement is available for the same subject under the same condition. Examples include repeated technical measurements, multiple regions or lesions within a patient, or several readings taken at the same visit. These measurements are replicate observations within a subject-condition cell. They should not be treated as independent subjects, because observations from the same subject are generally dependent. Replicate measurements are specified by the replicate argument. The replicate variable identifies repeated observations within the same subject-condition cell.

nparLD(response ~ group * time,
       data = dat,
       subject = "subject",
       replicate = "replicate")

For hypotheses in relative marginal effects, the argument cell.weights determines how replicate measurements contribute to the estimated effects.

With

cell.weights = "subjects"

replicate scores are first averaged within each observed subject-condition cell. The resulting effect describes a typical subject-condition cell. This is the natural choice when subjects are the primary sampling units and replicate measurements are repeated observations within those units.

With

cell.weights = "observations"

all replicate observations contribute directly. Subject-condition cells with more replicate observations therefore receive more weight. This targets a typical replicate observation and is appropriate only when the replicate observations themselves are part of the scientific target.

For hypotheses in marginal distribution functions, dependent replicates are handled through subject-level score averaging. Thus, inference remains based on the independent subject-level units.


data(brdu)

fit_brdu_subjects <- nparLD(
  resp ~ dose,
  data = brdu,
  subject = "culture",
  replicate = "replicate",
  hypothesis = "H0p",
  cell.weights = "subjects"
)


fit_brdu_observations <- nparLD(
  resp ~ dose,
  data = brdu,
  subject = "culture",
  replicate = "replicate",
  hypothesis = "H0p",
  cell.weights = "observations"
)

fit_brdu_subjects
#> 
#> Nonparametric Longitudinal Analysis
#> -------------------------------------------------
#> Design: LD-F1-Design 
#> 
#> Subplot factors: dose
#> Number of subjects: 5
#> Hypotheses: Tested in marginal Mann-Whitney effects (relative marginal effects)
#> Ranking: Pseudo-ranks
#> Confidence interval method: logit
#> 
#> Estimated relative effects
#> -------------------------------------------------
#>   dose Nsubj lambda Nobs Nmiss ScoreMean    RTE std.error  lower  upper
#> 1    0     5      5   15     0     22.13 0.3606   0.10454 0.2748 0.4562
#> 2  0.1     5      5   15     0     17.00 0.2750   0.06646 0.2207 0.3369
#> 3    1     5      5   15     0     38.70 0.6367   0.12584 0.5210 0.7384
#> 4   10     5      5   15     0     44.17 0.7278   0.17464 0.5525 0.8527
#> 
#> Global test results
#> -------------------------------------------------
#> 
#> Wald-type statistic (WTS)
#>      Statistic df p-value
#> dose  614.6248  3       0
#> 
#> ANOVA-type statistic (ATS)
#>      Statistic    df p-value
#> dose   11.4002 1.623   1e-04
fit_brdu_observations
#> 
#> Nonparametric Longitudinal Analysis
#> -------------------------------------------------
#> Design: LD-F1-Design 
#> 
#> Subplot factors: dose
#> Number of subjects: 5
#> Hypotheses: Tested in marginal Mann-Whitney effects (relative marginal effects)
#> Ranking: Pseudo-ranks
#> Confidence interval method: logit
#> 
#> Estimated relative effects
#> -------------------------------------------------
#>   dose Nsubj lambda Nobs Nmiss ScoreMean    RTE std.error  lower  upper
#> 1    0     5      5   15     0     22.13 0.3606   0.10454 0.2748 0.4562
#> 2  0.1     5      5   15     0     17.00 0.2750   0.06646 0.2207 0.3369
#> 3    1     5      5   15     0     38.70 0.6367   0.12584 0.5210 0.7384
#> 4   10     5      5   15     0     44.17 0.7278   0.17464 0.5525 0.8527
#> 
#> Global test results
#> -------------------------------------------------
#> 
#> Wald-type statistic (WTS)
#>      Statistic df p-value
#> dose  614.6248  3       0
#> 
#> ANOVA-type statistic (ATS)
#>      Statistic    df p-value
#> dose   11.4002 1.623   1e-04

For relative marginal effects, cell.weights = "subjects" targets a typical subject-condition cell, while cell.weights = "observations" targets a typical replicate observation. In the example, the numbers of replicates are identical and therefore, the results are identical.

Term-specific effect plots

The default plot method displays the estimated cell-level relative effects.

plot(fit_shoulder)

If factor-specific information is requested by setting Factor.Information = TRUE, the plot method can also display estimated relative effects and confidence intervals for selected model terms. This is useful when the scientific question concerns a main effect or interaction rather than all factorial cells.

fit_shoulder_info <- nparLD(
  resp ~ group1 * group2 * time,
  data = shoulder,
  subject = "subject",
  hypothesis = "H0p",
  Factor.Information = TRUE
)

plot(fit_shoulder_info, term = "group1:time")

Several terms can be displayed in one faceted plot.

plot(fit_shoulder_info, term = c("time", "group1:time"))

The horizontal reference line at 0.5 marks the value corresponding to no tendency relative to the reference distribution. Values above 0.5 indicate a tendency toward larger responses, whereas values below 0.5 indicate a tendency toward smaller responses. The plots are ggplot2 plots and hence can be adapted (e.g., in crease font) in the usual way.

Test statistics and multiple contrasts

For each specified hypothesis, nparLD() returns Wald-type and ANOVA-type test statistics. Both statistics are based on the estimated (relative marginal Mann-Whitney) effect vector and its estimated covariance matrix, but they use this information differently. Test procedures testing "H0F" estimate the covariance matrix under this null hypothesis, whereas all procedures testing "H0p" use a covariance matrix estimator that is consistent under arbitrary (but fixed) alternatives.

The Wald-type statistic uses the full covariance structure of the estimated effects. In general form, it can be written as

\[ Q_N = N \widehat{\mathbf{\theta}}^\top C^\top (\mathbf{C} \widehat{\mathbf{V}}_N\mathbf{C}^\top)^+ \mathbf{C} \widehat{\mathbf{\theta}}, \]

where \(\widehat{\mathbf{\theta}}\) denotes the estimated target vector, \(\widehat{\mathbf{V}}_N\) is the estimated covariance matrix, \(\mathbf{C}\) is the contrast matrix defining the hypothesis, and \((\cdot)^+\) denotes a generalized inverse.

The ANOVA-type statistic is based on a quadratic form using the projection matrix associated with the hypothesis,

\[ A_N = \frac{ N \widehat{\mathbf{\theta}}^\top T \widehat{\mathbf{\theta}} }{ \operatorname{tr}(\mathbf{T} \widehat{\mathbf{V}}_N)}. \]

It provides an alternative approximation that is often used in nonparametric factorial methods, especially when sample sizes are small.

In addition to global tests for main effects and interactions, nparLD() can compute multiple contrast procedures for selected factors or interactions. The argument contrast specifies the effect for which local comparisons should be computed. For a selected factor or interaction, let \(\mathbf{K}\) denote the contrast matrix used for the multiple contrast procedure. The estimated local contrasts are

\[ \widehat{\mathbf{\delta}} = \mathbf{K} \widehat{\mathbf{\theta}}, \]

with covariance estimator

\[ \widehat{\mathbf{\Sigma}}_{\delta} = \mathbf{K} \widehat{\mathbf{V}}_N \mathbf{K}^\top . \]

The corresponding studentized contrast statistics are

\[ T_\ell = \frac{ \sqrt{N}\widehat{\delta}_\ell }{ \sqrt{\widehat{\Sigma}_{\delta,\ell\ell}} }, \qquad \ell = 1,\ldots,q, \]

where \(q\) is the number of local contrasts. Simultaneous confidence intervals are obtained from the joint distribution of these studentized statistics, using the estimated correlation matrix derived from \(\widehat{\mathbf{\Sigma}}_{\delta}\).

For example,

contrast = list("time", "Dunnett")

requests Dunnett-type contrasts for the levels of time, while

contrast = list("group:time")

requests a multiple contrast procedure for the interaction between group and time.

When multiple contrasts are requested, the output contains a global multiple contrast test and local contrast estimates with standard errors, simultaneous confidence intervals, test statistics, adjusted p-values, and degrees of freedom. The contrast results can also be plotted:

plot(fit$MCTP)

Note that all procedures testing the null hypothesis "H0F" can not yield confidence intervals for the relative marginal effects.

A real-data example: shoulder tip pain

The shoulder data illustrate a factorial longitudinal design with two whole-plot factors and one repeated-measures factor. The response variable resp is a shoulder pain score measured repeatedly over time. The variables group1 and group2 define between-subject factors, while time is the within-subject factor.

data(shoulder)

str(shoulder)
#> 'data.frame':    246 obs. of  5 variables:
#>  $ resp   : int  1 3 1 1 2 1 3 2 1 2 ...
#>  $ time   : int  1 1 1 1 1 1 1 1 1 1 ...
#>  $ group1 : Factor w/ 2 levels "N","Y": 2 2 2 2 2 2 2 2 2 2 ...
#>  $ group2 : Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
#>  $ subject: int  1 3 4 5 8 9 10 12 16 18 ...

This corresponds to an F2-LD-F1 design in the classical terminology: two whole-plot factors and one longitudinal factor. In the redesigned interface, the model is specified directly by the formula

fit_shoulder_F <- nparLD(
  resp ~ group1 * group2 * time,
  data = shoulder,
  subject = "subject",
  hypothesis = "H0F"
)

fit_shoulder_F
#> 
#> Nonparametric Longitudinal Analysis
#> -------------------------------------------------
#> Design: F2-LD-F1-Design 
#> 
#> Whole-plot factors: group1, group2
#> Subplot factors: time
#> Number of subjects: 41
#> Hypotheses: Tested in distribution functions
#> Ranking: Pseudo-ranks
#> 
#> Estimated relative effects
#> -------------------------------------------------
#>    group1 group2 time Nsubj lambda Nobs Nmiss ScoreMean    RTE std.error
#> 1       N      F    1    11     11   11     0    152.71 0.6187    0.6079
#> 2       N      F    2    11     11   11     0    172.57 0.6995    0.5054
#> 3       N      F    3    11     11   11     0    160.56 0.6507    0.5208
#> 4       N      F    4    11     11   11     0    180.86 0.7332    0.4936
#> 5       N      F    5    11     11   11     0    144.72 0.5863    0.5512
#> 6       N      F    6    11     11   11     0    130.89 0.5301    0.4501
#> 7       N      M    1     8      8    8     0    124.69 0.5048    0.6199
#> 8       N      M    2     8      8    8     0    166.79 0.6760    0.6101
#> 9       N      M    3     8      8    8     0    175.24 0.7103    0.6517
#> 10      N      M    4     8      8    8     0    171.01 0.6931    0.6326
#> 11      N      M    5     8      8    8     0    148.38 0.6011    0.6784
#> 12      N      M    6     8      8    8     0    120.46 0.4877    0.5643
#> 13      Y      F    1    14     14   14     0    121.82 0.4932    0.4429
#> 14      Y      F    2    14     14   14     0     98.47 0.3983    0.4139
#> 15      Y      F    3    14     14   14     0     87.35 0.3530    0.3378
#> 16      Y      F    4    14     14   14     0     99.22 0.4013    0.3602
#> 17      Y      F    5    14     14   14     0     70.29 0.2837    0.1545
#> 18      Y      F    6    14     14   14     0     82.16 0.3319    0.2461
#> 19      Y      M    1     8      8    8     0    105.22 0.4257    0.5306
#> 20      Y      M    2     8      8    8     0    110.75 0.4482    0.4712
#> 21      Y      M    3     8      8    8     0     85.13 0.3440    0.3540
#> 22      Y      M    4     8      8    8     0     74.74 0.3018    0.2703
#> 23      Y      M    5     8      8    8     0     89.98 0.3637    0.4470
#> 24      Y      M    6     8      8    8     0     89.98 0.3637    0.4470
#> 
#> Global test results
#> -------------------------------------------------
#> 
#> Wald-type statistic (WTS)
#>                    Statistic df p-value
#> group1               16.5421  1  0.0000
#> group2                0.0474  1  0.8277
#> time                 17.3493  5  0.0039
#> group1:group2         0.0320  1  0.8581
#> group1:time          28.9628  5  0.0000
#> group2:time          12.5250  5  0.0283
#> group1:group2:time    5.1069  5  0.4030
#> 
#> ANOVA-type statistic (ATS)
#>                    Statistic     df p-value
#> group1               16.5421 1.0000  0.0000
#> group2                0.0474 1.0000  0.8277
#> time                  3.4781 2.7007  0.0188
#> group1:group2         0.0320 1.0000  0.8581
#> group1:time           3.8134 2.7007  0.0123
#> group2:time           1.1656 2.7007  0.3194
#> group1:group2:time    0.4371 2.7007  0.7057

The analysis with hypothesis = "H0F" tests hypotheses in the marginal distribution functions. These tests compare the complete marginal distributions defined by the factorial cells and may therefore detect differences beyond location shifts.

For effect-based interpretation, we can formulate the analysis in terms of unweighted relative marginal effects:

fit_shoulder_p <- nparLD(
  resp ~ group1 * group2 * time,
  data = shoulder,
  subject = "subject",
  hypothesis = "H0p"
)

fit_shoulder_p
#> 
#> Nonparametric Longitudinal Analysis
#> -------------------------------------------------
#> Design: F2-LD-F1-Design 
#> 
#> Whole-plot factors: group1, group2
#> Subplot factors: time
#> Number of subjects: 41
#> Hypotheses: Tested in marginal Mann-Whitney effects (relative marginal effects)
#> Ranking: Pseudo-ranks
#> Confidence interval method: logit
#> 
#> Estimated relative effects
#> -------------------------------------------------
#>    group1 group2 time Nsubj lambda Nobs Nmiss ScoreMean    RTE std.error  lower
#> 1       N      F    1    11     11   11     0    152.71 0.6187    0.5560 0.4410
#> 2       N      F    2    11     11   11     0    172.57 0.6995    0.4518 0.5466
#> 3       N      F    3    11     11   11     0    160.56 0.6507    0.4547 0.5024
#> 4       N      F    4    11     11   11     0    180.86 0.7332    0.4278 0.5845
#> 5       N      F    5    11     11   11     0    144.72 0.5863    0.5004 0.4297
#> 6       N      F    6    11     11   11     0    130.89 0.5301    0.4050 0.4068
#> 7       N      M    1     8      8    8     0    124.69 0.5048    0.5573 0.3400
#> 8       N      M    2     8      8    8     0    166.79 0.6760    0.5067 0.5068
#> 9       N      M    3     8      8    8     0    175.24 0.7103    0.5629 0.5149
#> 10      N      M    4     8      8    8     0    171.01 0.6931    0.5320 0.5123
#> 11      N      M    5     8      8    8     0    148.38 0.6011    0.5736 0.4202
#> 12      N      M    6     8      8    8     0    120.46 0.4877    0.4606 0.3512
#> 13      Y      F    1    14     14   14     0    121.82 0.4932    0.4350 0.3636
#> 14      Y      F    2    14     14   14     0     98.47 0.3983    0.3976 0.2848
#> 15      Y      F    3    14     14   14     0     87.35 0.3530    0.3388 0.2573
#> 16      Y      F    4    14     14   14     0     99.22 0.4013    0.3633 0.2967
#> 17      Y      F    5    14     14   14     0     70.29 0.2837    0.2157 0.2225
#> 18      Y      F    6    14     14   14     0     82.16 0.3319    0.2970 0.2480
#> 19      Y      M    1     8      8    8     0    105.22 0.4257    0.5490 0.2715
#> 20      Y      M    2     8      8    8     0    110.75 0.4482    0.5027 0.3036
#> 21      Y      M    3     8      8    8     0     85.13 0.3440    0.3695 0.2411
#> 22      Y      M    4     8      8    8     0     74.74 0.3018    0.2831 0.2227
#> 23      Y      M    5     8      8    8     0     89.98 0.3637    0.4337 0.2437
#> 24      Y      M    6     8      8    8     0     89.98 0.3637    0.4337 0.2437
#>     upper
#> 1  0.7695
#> 2  0.8180
#> 3  0.7746
#> 4  0.8429
#> 5  0.7271
#> 6  0.6498
#> 7  0.6686
#> 8  0.8090
#> 9  0.8500
#> 10 0.8293
#> 11 0.7581
#> 12 0.6260
#> 13 0.6237
#> 14 0.5238
#> 15 0.4622
#> 16 0.5157
#> 17 0.3541
#> 18 0.4281
#> 19 0.5958
#> 20 0.6021
#> 21 0.4640
#> 22 0.3947
#> 23 0.5036
#> 24 0.5036
#> 
#> Global test results
#> -------------------------------------------------
#> 
#> Wald-type statistic (WTS)
#>                    Statistic df p-value
#> group1               16.9539  1  0.0000
#> group2                0.0476  1  0.8274
#> time                 14.7195  5  0.0116
#> group1:group2         0.0318  1  0.8584
#> group1:time          22.6602  5  0.0004
#> group2:time          12.1939  5  0.0322
#> group1:group2:time    5.2195  5  0.3897
#> 
#> ANOVA-type statistic (ATS)
#>                    Statistic     df p-value
#> group1               16.9539 1.0000  0.0000
#> group2                0.0476 1.0000  0.8274
#> time                  3.3876 2.7526  0.0203
#> group1:group2         0.0318 1.0000  0.8584
#> group1:time           3.7277 2.7325  0.0133
#> group2:time           1.1548 2.7092  0.3234
#> group1:group2:time    0.4311 2.6790  0.7084

The estimated relative effects can be displayed graphically:

plot(fit_shoulder_p)

A value above 0.5 indicates a tendency toward larger responses relative to the unweighted reference distribution, while a value below 0.5 indicates a tendency toward smaller responses.

Finally, we can request multiple contrasts for an interaction. For example, the following call computes contrasts for the interaction between group1 and time:

fit_shoulder_contrast <- nparLD(
  resp ~ group1 * group2 * time,
  data = shoulder,
  subject = "subject",
  hypothesis = "H0p",
  contrast = list("group1:time")
)

fit_shoulder_contrast$MCTP
#> 
#> Multiple contrast test procedure
#> -------------------------------------------------
#> Factor or interaction: group1:time
#> SCI method: fisher
#> 
#> Global result
#> -------------------------------------------------
#>   Statistic p.value
#> 1    3.5025  0.0103
#> 
#> Local contrast results
#> -------------------------------------------------
#>      Estimate Std.Error   lower   upper Statistic p.value df
#> C 1   -0.0731    0.2123 -0.1640  0.0190   -2.2089  0.1552 21
#> C 2    0.0080    0.1538 -0.0588  0.0747    0.3317  0.9970 21
#> C 3    0.0417    0.1233 -0.0119  0.0950    2.1667  0.1675 21
#> C 4    0.0565    0.1034  0.0116  0.1012    3.5025  0.0103 21
#> C 5    0.0107    0.1436 -0.0517  0.0730    0.4772  0.9850 21
#> C 6   -0.0438    0.1467 -0.1071  0.0199   -1.9124  0.2593 21
#> C 7    0.0731    0.2123 -0.0190  0.1640    2.2089  0.1553 21
#> C 8   -0.0080    0.1538 -0.0747  0.0588   -0.3317  0.9970 21
#> C 9   -0.0417    0.1233 -0.0950  0.0119   -2.1667  0.1674 21
#> C 10  -0.0565    0.1034 -0.1012 -0.0116   -3.5025  0.0105 21
#> C 11  -0.0107    0.1436 -0.0730  0.0517   -0.4772  0.9850 21
#> C 12   0.0438    0.1467 -0.0199  0.1071    1.9124  0.2591 21
#> 
#> Contrast matrix: available in component `Contrast.Matrix`.
#> Use `print(x, show.matrix = TRUE)` to display it.
plot(fit_shoulder_contrast$MCTP)

This output contains both a global multiple contrast test and local contrast results with simultaneous confidence intervals.

Interpreting the output

The object returned by nparLD() contains several components. Printing the object gives a compact summary of the design, the selected hypothesis type, the estimated effects, and the available test results. The most commonly used components are effects, WTS, ATS, and, when multiple contrasts are requested, MCTP.

The component effects contains the estimated relative effects for the factorial cells. It includes the factor-level combinations, the number of contributing subjects and observations, the number of missing observations, the estimated relative effect, and its standard error. For analyses with hypothesis = "H0p", confidence limits are also returned.

head(fit_shoulder_p$effects)
#>   group1 group2 time Nsubj lambda Nobs Nmiss ScoreMean       RTE std.error
#> 1      N      F    1    11     11   11     0  152.7108 0.6187432 0.5559869
#> 2      N      F    2    11     11   11     0  172.5747 0.6994909 0.4517811
#> 3      N      F    3    11     11   11     0  160.5610 0.6506543 0.4546888
#> 4      N      F    4    11     11   11     0  180.8620 0.7331790 0.4278484
#> 5      N      F    5    11     11   11     0  144.7200 0.5862603 0.5003957
#> 6      N      F    6    11     11   11     0  130.8940 0.5300571 0.4050453
#>       lower     upper
#> 1 0.4409740 0.7695275
#> 2 0.5466137 0.8179845
#> 3 0.5024039 0.7745558
#> 4 0.5845221 0.8429388
#> 5 0.4297327 0.7271074
#> 6 0.4067651 0.6497866

The components WTS and ATS contain the global tests for the main effects and interactions defined by the model formula.

fit_shoulder_p$WTS
#>                      Statistic df      p-value
#> group1             16.95386542  1 3.829926e-05
#> group2              0.04755604  1 8.273719e-01
#> time               14.71947498  5 1.163070e-02
#> group1:group2       0.03183416  1 8.583921e-01
#> group1:time        22.66020930  5 3.919732e-04
#> group2:time        12.19393660  5 3.222489e-02
#> group1:group2:time  5.21951962  5 3.896815e-01
fit_shoulder_p$ATS
#>                      Statistic       df      p-value
#> group1             16.95386542 1.000000 3.829926e-05
#> group2              0.04755604 1.000000 8.273719e-01
#> time                3.38759671 2.752603 2.033849e-02
#> group1:group2       0.03183416 1.000000 8.583921e-01
#> group1:time         3.72770726 2.732513 1.334084e-02
#> group2:time         1.15480930 2.709224 3.234477e-01
#> group1:group2:time  0.43105480 2.678980 7.083611e-01

The Wald-type statistic uses the estimated covariance matrix directly, whereas the ANOVA-type statistic uses a trace-based approximation. In practice, both statistics are useful summaries of the global hypotheses. Differences between the two may occur, especially in small samples, unbalanced designs, or settings with singular covariance structures.

When multiple contrasts are requested, the component MCTP contains the corresponding contrast results.

fit_shoulder_contrast$MCTP
#> 
#> Multiple contrast test procedure
#> -------------------------------------------------
#> Factor or interaction: group1:time
#> SCI method: fisher
#> 
#> Global result
#> -------------------------------------------------
#>   Statistic p.value
#> 1    3.5025  0.0103
#> 
#> Local contrast results
#> -------------------------------------------------
#>      Estimate Std.Error   lower   upper Statistic p.value df
#> C 1   -0.0731    0.2123 -0.1640  0.0190   -2.2089  0.1552 21
#> C 2    0.0080    0.1538 -0.0588  0.0747    0.3317  0.9970 21
#> C 3    0.0417    0.1233 -0.0119  0.0950    2.1667  0.1675 21
#> C 4    0.0565    0.1034  0.0116  0.1012    3.5025  0.0103 21
#> C 5    0.0107    0.1436 -0.0517  0.0730    0.4772  0.9850 21
#> C 6   -0.0438    0.1467 -0.1071  0.0199   -1.9124  0.2593 21
#> C 7    0.0731    0.2123 -0.0190  0.1640    2.2089  0.1553 21
#> C 8   -0.0080    0.1538 -0.0747  0.0588   -0.3317  0.9970 21
#> C 9   -0.0417    0.1233 -0.0950  0.0119   -2.1667  0.1674 21
#> C 10  -0.0565    0.1034 -0.1012 -0.0116   -3.5025  0.0105 21
#> C 11  -0.0107    0.1436 -0.0730  0.0517   -0.4772  0.9850 21
#> C 12   0.0438    0.1467 -0.0199  0.1071    1.9124  0.2591 21
#> 
#> Contrast matrix: available in component `Contrast.Matrix`.
#> Use `print(x, show.matrix = TRUE)` to display it.

The local contrast table contains the estimated contrasts, standard errors, simultaneous confidence limits, test statistics, adjusted p-values, and degrees of freedom. These results are most naturally interpreted on the relative-effect scale. A simultaneous confidence interval that does not contain zero indicates a contrast that is significant after adjustment for the selected family of comparisons.

The estimated covariance matrix can be included in the output by setting covariance = TRUE. This option is mainly intended for methodological checks, advanced post-processing, or reproducibility of the reported test statistics.

For graphical summaries, plot() can be applied to the fitted object or to the multiple contrast component.

plot(fit_shoulder_p)
plot(fit_shoulder_contrast$MCTP)

Practical recommendations

The redesigned nparLD() function provides several analysis options. The following recommendations may help to choose an appropriate specification.

Use hypothesis = "H0F" when the primary question concerns equality of marginal distribution functions. These hypotheses compare complete marginal distributions. This option is useful when the scientific question is distributional and not restricted to a single effect measure. Classical rank procedures test null hypotheses in distribution functions.

Use hypothesis = "H0p" when the primary question concerns relative marginal effects. These hypotheses are formulated for unweighted relative effects and are therefore particularly useful for effect interpretation, graphical summaries, multiple contrasts, and simultaneous confidence intervals. The resulting estimates can be interpreted on the relative-effect scale: values above 0.5 indicate a tendency toward larger responses than the reference distribution, whereas values below 0.5 indicate a tendency toward smaller responses.

Weighted relative effects should be interpreted as descriptive summaries. They are based on the sample-size weighted reference distribution and are closely connected to ordinary ranks. Because they depend on the observed sample-size allocation, they are not used as the inferential target for hypotheses in relative marginal effects.

Unweighted relative effects should be used when inference on relative effects is desired. They are based on an equally weighted reference distribution over the factorial cells and are estimated by pseudo-ranks. This keeps the interpretation of the effects separate from the sample-size allocation, which is especially important in unbalanced designs.

For incomplete longitudinal data, missing response values should be coded as NA, while the subject identifier and the relevant factor-level information should remain in the data set. This allows nparLD() to identify the incomplete subject-condition cells automatically. Subjects with partially observed trajectories are not discarded completely as long as they contribute observations to at least part of the design.

For dependent replicate measurements, use the replicate argument to identify multiple observations within the same subject-condition cell. With cell.weights = "subjects", replicate scores are averaged within subject-condition cells, so that the analysis targets a typical subject-condition cell. With cell.weights = "observations", all replicate observations contribute directly, so that the analysis targets a typical replicate observation.

In most applications, a useful workflow is to first identify the primary research questions of interest and translate them into suitable contrasts. The multiple contrast procedures can then be used for selected factors or interactions to obtain local comparisons and simultaneous confidence intervals. Graphical displays of the estimated effects and simultaneous confidence intervals can help to communicate the results.

Global tests in WTS and ATS are useful when the global null hypothesis is the primary or sole object of interest. In contrast, when the scientific question concerns specific group, time, or interaction comparisons, contrast procedures usually provide the more directly interpretable results.