
---
title: "Nonparametric Longitudinal Data Analysis with nparLD"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Nonparametric Longitudinal Data Analysis with nparLD}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```


```{r setup}
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


```{r, eval=FALSE}

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

```{r, eval = FALSE}
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,

```{r, eval = FALSE}
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

```{r dental-example}
data(dental)

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

fit_dental
```

## A factorial longitudinal example

```{r shoulder-example, fig.width = 7, fig.height = 5}
data(shoulder)

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

fit_shoulder
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.


```{r missing-example}
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
```
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.

```{r, eval = FALSE}
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

```{r, eval = FALSE}
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

```{r, eval = FALSE}
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.


```{r replicate-example}

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
fit_brdu_observations

```

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.

```{r plot-cell-effects, fig.width=7, fig.height=4}
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.


```{r plot-term effects, fig.width=7, fig.height=4}
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.

```{r plot-cell-effectsseveral, fig.width=7, fig.height=4}
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,

```{r, eval = FALSE}
contrast = list("time", "Dunnett")
```

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

```{r, eval = FALSE}
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:

```{r, eval = FALSE}
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.

```{r shoulder-data}
data(shoulder)

str(shoulder)
```

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

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

fit_shoulder_F
```

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:

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

fit_shoulder_p
```

The estimated relative effects can be displayed graphically:

```{r shoulder-plot, fig.width = 7, fig.height = 5}
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`:

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

fit_shoulder_contrast$MCTP
```

```{r shoulder-mctp-plot, fig.width = 7, fig.height = 5}
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.

```{r effects-component}
head(fit_shoulder_p$effects)
```

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

```{r test-components}
fit_shoulder_p$WTS
fit_shoulder_p$ATS
```

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.

```{r mctp-component}
fit_shoulder_contrast$MCTP
```

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.

```{r output-plots, eval = FALSE}
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.




