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