sanon : An R Package for Stratified Analysis with Nonparametric Covariable Adjustment

Kawaguchi, Koch, and Wang (2011) provide methodology and applications for a strat-iﬁed Mann-Whitney estimator that addresses the same comparison between two randomized groups for a strictly ordinal response variable as the van Elteren test statistic for randomized clinical trials with strata. The sanon package provides the implementation of the method within the R programming environment. The usage of sanon is illustrated with ﬁve examples. The ﬁrst example is a randomized clinical trial with eight strata and a univariate ordinal response variable. The second example is a randomized clinical trial with four strata, two covariables, and four ordinal response variables. The third example is a crossover design randomized clinical trial with two strata, one covariable, and two ordinal response variables. The fourth example is a randomized clinical trial with seven strata (which are managed as a categorical covariable), three ordinal covariables with missing values, and three ordinal response variables with missing values. The ﬁfth example is a randomized clinical trial with six strata, a categorical covariable with three levels, and three ordinal response variables with missing values.


Introduction
The primary analyses for confirmatory randomized clinical trials (and particularly those with regulatory objectives) should consist of protocol specified methods that have minimal assumptions.A nonparametric approach such as the Wilcoxon (or Mann-Whitney) test statistic for the comparison between two treatments through the ranking of a response variable for all patients (in the pooled treatment groups) has essentially no assumptions (under the strong null hypothesis of no treatment differences in the sense that each patient has the same response to both treatments) beyond valid randomization as its basis.More generally, Mann-Whitney test statistics and the corresponding estimators address the global null hypothesis of equality of response distributions for two treatment groups through their implied null hypothesis of P(Group 1 response > Group 2 response) = ξ = 0.5; and their sensitivity and power to detect differences between two distributions depend jointly on the extent to which |ξ − 0.5| > 0 and the applicable sample size.Nevertheless, such methods have the recognized limitation that their power for comparing different distributions for which ξ = 0.5 equals the specified type 1 error regardless of the applicable sample sizes (although differences between such distributions are typically of relatively minimal interest for clinical trials that evaluate whether one treatment group has relatively more patients with better responses than the other treatment group).The coin package (Hothorn, Hornik, Wiel, and Zeileis 2006, 2008, 2015) in the R programming environment (R Core Team 2015) can handle ordered and multivariate responses as well as stratification and provides both randomization/permutation tests and asymptotic tests based on conditional inference.NParCov3 (Zink and Koch 2012) is a SAS/IML macro (SAS Institute Inc. 2011) written to conduct the nonparametric randomization-based covariance analyses of Koch, Tangen, Jung, and Amara (1998).
For confirmatory randomized clinical trials with stratified designs for comparing two treatments, Kawaguchi et al. (2011) propose stratified multivariate Mann-Whitney estimators as a useful structure for the analysis of strictly ordinal response variables.Their scope can address strata with at least minimal sample sizes (e.g., ≥ 16), and randomization-based covariance adjustment is possible.The method is based on the Mann-Whitney estimator for the probability that a randomly selected patient from one treatment group has better status for a response variable within a stratum than a randomly selected patient from the other treatment group (with ties being randomly broken with probability 0.5).Such Mann-Whitney estimators can be combined across the strata to provide stratified estimator counterparts that address the same comparisons between the two treatment groups as the van Elteren test statistic.The multivariate vector of such stratified Mann-Whitney estimators for multivariate response variables can be considered for one or more response variables such as in repeated measurements and these can have missing completely at random (MCAR) data (or missing data for which direct imputation methods are applicable).Randomization-based covariable adjustment is possible for stratified multivariate Mann-Whitney estimators by expanding the vector of such estimators to include stratified differences between means of covariables.The latter estimators for the covariables then have constraints to 0's invoked.For this purpose, weighted least squares methods are applied with weights based on the estimated covariance matrix for the expanded vector from the methods for ratios of multivariate U-statistics; see Stokes, Davis, and Koch (2012, Chapter 14) and Koch et al. (1998).The resulting estimators are stratified multivariate Mann-Whitney estimators with randomization-based covariance adjustment, and their interpretation is the same as the original Mann-Whitney estimator mentioned above.With sufficiently large sample sizes, such estimators have an approximately multivariate normal distribution with the covariance matrix being essentially known through its corresponding consistent estimator.Accordingly, confidence intervals can be constructed for linear functions of such fully adjusted Mann-Whitney estimators (with respect to both stratification and covariables), and the scope of such linear functions can include the separate response variables, averages across response variables, and contrasts among response variables.
The software implementing the method was developed in the R programming environment.The sanon package (Kawaguchi 2015) contains functions to apply the method as well as two example data sets which were used in Kawaguchi et al. (2011) and three data sets from other papers.The main function in the package has a similar structure as functions in R for standard statistical methods such as lm for linear models.This paper illustrates the usage of the package with five example data sets.Section 2 describes the methods, Section 3 explains the code, and Section 5 illustrates those methods for five examples introduced in Section 4.

Methods
As previously outlined in Section 1, a stratified Mann-Whitney estimator and corresponding confidence interval can be constructed to address the same comparison between two randomized groups as the van Elteren test statistic.The specifications for the formal structure for this method are in Section 2.1 for a multivariate set of r response variables.Section 2.2 discusses randomization-based covariance adjustment for stratified multivariate Mann-Whitney estimators and explains the corresponding constraints for no expected differences between randomized groups for stratified means of covariables.

Stratified multivariate Mann-Whitney estimator
Let h = 1, 2, . . ., q index a set of strata within which patients are randomized to two groups indexed by i = 1, 2. Let k = 1, 2, . . ., r index the response variables with observations for the n hi patients in the ith group of the hth stratum; and some of these response variables can be baselines at times prior to any treatment for a group; or in crossover studies, they can be at washout times prior to the treatments which are subsequent to the first treatment.Let j = 1, 2, . . ., N index the patients after pooling of all patients in the clinical trial regardless of their groups or strata.Let S j denote the stratum for the jth patient, and let t j correspond to the group for the jth patient with t j = 1 if i = 1 for patient j and t j = −1 if i = 2 for patient j.Let Y j = (Y j1 , . . ., Y jr ) denote the response vector for the jth patient with Y jk denoting the kth strictly ordinal response variable for the jth patient.Since some of the Y jk may be missing (by a MCAR process), let Z jk = 1 if Y jk is not missing and let Z jk = 0 if Y jk is missing; and let Z j = (Z j1 , . . ., Z jr ) .In the subsequent discussion, any missing Y jk operationally has a replacement by 0 since the value used for such replacement has no role in the subsequently described processes for estimation.
The stratified Mann-Whitney estimator for the kth response variable is ξk = ( θ1k / θ2k ) for which θ1k and θ2k are defined in (1) and ( 2 (2) Also, I(A) has the value 1 if the condition A is satisfied or the value 0 otherwise, and N jk denotes the sample size for the kth response variable for patients with the same stratum and group as the jth patient.For (1) and ( 2), the numerator of U 1jj k equals 1 for a pair of patients who are in the same stratum with different groups and for whom the kth response variable is larger for the patient in group 1 than for the patient in group 2, and it equals 0 for all other pairs of patients; and the numerator of U 2jj k equals 1 for a pair of patients who are in the same stratum with different groups and for whom the kth response variable is observed (i.e., not missing), and it equals 0 for all other pairs of patients.Accordingly, the numerator of θ2k equals q h=1 w hk for which w hk = (n h1k n h2k )/(n h * k + 1) with n hik being the number of patients for whom the kth response variable is observed among the n hi patients for the ith group in the hth stratum and n h * k = (n h1k + n h2k ); and the numerator of θ1k equals q h=1 w hk ξhk for which ξhk for the respective strata are the within strata Mann-Whitney estimators for the proportions of pairs of patients with different groups and for whom the kth response variable is larger for the group 1 patient than the group 2 patient.Thus, ξk = q h=1 w hk ξhk q h=1 w hk is a weighted average of the within strata, Mann-Whitney estimators ξhk with the w hk as weights.
) denote a compound vector for the jth patient where U 1j = (U 1j1 , . . ., U 1jr ) and U 2j = (U 2j1 , . . ., U 2jr ) .Let F = N j=1 F j /N denote the sample mean vector for the F j .As noted in Davis and Quade (1968), Puri and Sen (1971), Quade (1974), Carr, Hafner, andKoch (1989), Jung andKoch (1998), andJung andKoch (1999), a consistent estimator for the covariance matrix for F is given in (3). (3) The stratified Mann-Whitney estimator for ξ = (ξ 1 , ξ 2 , . . ., ξ r ) is ξ = D −1 θ2 θ1 = ( ξ1 , ξ2 , . . ., ξr ) , where D a denotes a diagonal matrix with the elements of vector a on the main diagonal.Here θ1 = ( θ11 , θ12 , . . ., θ1r ) corresponds to the first r elements of F (as they represent numerators of the ξk for the r response variables); and θ2 = ( θ21 , θ22 , . . ., θ2r ) corresponds to the remaining r elements of F ; i.e., F = ( θ 1 , θ 2 ) .Based on the Taylor series linearization, a consistent estimator V ξ for the covariance matrix for ξ is given in (4). (4) As noted in Kawaguchi et al. (2011), with sufficient sample size for ξ to have an approximately multivariate normal distribution (e.g., n +ik ≥ 50 and all n hik ≥ 4), then a two-sided 100(1 − 2α)% confidence interval for the linear statistic c ξ for the comparison between the two groups is c ξ ± z α c V ξ c where z α is the 100(1 − α) percentile of the standard normal distribution with mean 0 and variance 1.In this regard, the interpretation of the results from assessments of such contrasts may need some caution since Mann-Whitney estimators do not have explicitly transitive relationships with one another, although approximately transitive relationships can sometimes be applicable to their corresponding log odds ψk = log e [ ξk /(1 − ξk )] through Bradley-Terry specifications for paired comparisons; see Jung and Koch (1999) and Kawaguchi and Koch (2010).

Randomization-based covariance adjustment
Let x j = (x j1 , . . ., x jM ) denote the vector of M numeric covariables for the jth patient with x jm denoting the numeric value of the mth covariable for the jth patient.All of the M covariables have observations prior to the randomization of the patients to the two groups.Also, any categorical covariable has been expressed as a set of indicator variables that correspond to all except one of its categories with the excluded category serving as a reference category; and any ordinal covariable is managed in the same way as the response variables via (1)-(3).All of the covariables are assumed to have no missing data as is often realistic (although it is possible to manage a covariable with substantial missing data as a categorical covariable and to manage minimal missing data for other covariables in the same way as for response variables).
The stratification adjusted mean difference estimator for the mth covariable is where w h = n h1 n h2 /(n h1 + n h2 ), and φ1m and φ2 are defined as in (6).
In ( 6), the U 1jj m and U 2jj have the structure shown in ( 7) where N j = n hi if patient j is from the hth stratum and is in the ith group.Also, the numerator of U 1jj m (for a pair of patients who are from the same stratum with different groups) equals the difference between the patient in group 1 and the patient in group 2 for the mth covariable, and it equals 0 for all other pairs of patients; and the numerator of U 2jj has a comparable definition with respect to the covariables x j as U 2jj k has for the kth response variable.
Let f = ( ξ , g ) where ξ = D −1 θ2 θ1 and g = φ1 / φ2 ; here φ1 = ( φ11 , . . ., φ1M ) .The first r elements of f address comparisons between the two groups for the r response variables with stratification adjusted Mann-Whitney estimators, and the last M elements address comparisons between the two groups for the M covariables with stratification adjusted differences between their corresponding means.Since g would be expected to be null on the basis of randomization of patients to the two treatment groups, randomization-based covariance adjustment of ξ is possible by fitting the model P = [I r , 0 rM ] to f by weighted least squares; see Koch et al. (1998) andLaVange, Durham, andKoch (2005).The weights are based on a consistent estimator for the covariance matrix of f which is derived as follows. Let j=1 G j /N denote the sample mean vector for the G j .A consistent estimator for the covariance matrix of G is given in (8).
A consistent estimator for the covariance matrix of f is V f = HV G H where H as shown in ( 9) is from Taylor series linearization.
where V ξg corresponds to the covariances of ξ with g and V g corresponds to the covariance matrix of g.A consistent estimator for the covariance matrix of b is V b in (11).
Additional models that address the variation of the elements of b across the r response variables can be fit by weighted least squares methods.One such model can invoke randomizationbased covariance adjustment for any strictly ordinal covariables among the response variables by having a structure like P with rows of 0's corresponding to the strictly ordinal covariables and rows of an identity matrix corresponding to the strictly ordinal response variables.Through the resulting estimator b adj with full adjustment for all covariables and the strata, confidence intervals can be determined for linear statistics c b adj .Similarly, test statistics are applicable to assessments of homogeneity across the response variables for the adjusted estimators b adj for the differences between the two groups.
Since the preceding discussion indicates different roles for strata and covariables, some clarifying comments for the distinction between them is useful.In this regard, strata typically correspond to the cross-classification of categorical factors for stratified randomization.However, sometimes only one or two such factors that more clearly need balanced allocation of treatments and better enable reduction of within stratum variance are pre-specified as the structure for strata in situations where the use of all such factors would produce too many strata with possibly uninformatively small sample sizes for rankings of responses within them (e.g., (n h1 + n h2 ) ≤ 15).In these cases, the other factors for stratification can be pre-specified as covariables.In contrast, covariables are typically pre-specified baseline variables which potentially have strong associations with response variables (even though they did not serve as a factor for stratification).For such baseline variables, randomization-based covariance adjustment via (10) eliminates the spurious influence of random imbalances of their means for the treatment groups, and it provides variance reduction for treatment comparisons via (11).Also, as many as ten baseline variables (with appropriate pre-specified clinical justification) can have randomization-based covariance adjustment via (10) for clinical trials with at least a moderately large, overall sample size (e.g., N ≥ 300).In some situations, a categorical baseline variable that is not a factor for stratified randomization has pre-specified adjustment through stratification in order to remove the influence of random imbalances between treatment groups for its distribution and to have variance reduction for treatment comparisons more directly than provided by covariance adjustment.In summary, baseline variables can be either factors for stratification or covariables, and so pre-specification of which applies is very important to avoid potential bias from their roles having an arbitrary nature.Aspects of the scope of such analyses and the corresponding models have more specific discussion through the examples in Section 5.

Management of missing data
In addition to the method described in Section 2.1 to maintain missing values, three other options are newly introduced in this paper.The first option invokes the last observation carried forward (LOCF) convention for the kernels of U-statistics in (2).For this method, the denominators in ( 2) are revised to (N j + N j + 1), and missing values for the numerators are replaced by the last observed value of that variable; that is, each U 1jj k for the nearest observed preceding k has its numerator carried forward to the missing k and ) with U 2jj as in ( 7).If there is a missing value in the first measurement, then this measurement is managed as tied with the (observed or missing) other measurements in the corresponding stratum; and so Also, this management is more generally applied to all measurements prior to the first observed measurement by carrying forward U 1jj 1 and U 2jj 1 .The second option is the LOCF method based on the observed value of Y .As noted in the first LOCF method, the missing first value or values prior to the first observed value are managed as tied with the (observed or missing) other measurements in the corresponding stratum.The third option for missing data manages missing values as tied with all other values in the same stratum.For this method, U 2jj k = (N j + N j ) U 2jj /(N j + N j + 1) with U 2jj as in ( 7) and U 1jj k is modified to (12).
Additionally, the method described in Sections 2.1 and 2.2 is applicable with the method of multiple imputation for all missing data, but methods for multiple imputation are beyond the scope of this paper, particularly for strictly ordinal response variables.Nevertheless, adjusted Mann-Whitney estimates (10) and corresponding covariance matrices (11) from the respective invocations of multiple imputation can be combined to produce average estimators and their estimated covariance matrix by the methods of Rubin ( 2004) (as available, for example, in the MIANALYZE procedure of SAS).

Code
The R package sanon is available from the Comprehensive R Archive Network (CRAN) at http://CRAN.R-project.org/package=sanonwhere it can be installed from using the command install.packages("sanon").Then, you need to type library("sanon") to load the package.The main function sanon can be used by suitably specifying the input and it is formula based.The usage is: sanon(formula, data, P = NULL, res.na.action = "default") Argument Description formula A formula object, with the response on the left of a ~operator, and the terms on the right, in which functions strt, grp, covar, and catecovar are required to specify the role of variables on the right hand side.data A data.frame in which to interpret the variables named in the formula.
This function has four types of inputs (arguments), which are described in Table 1 and the details are below.
The formula consists of variable names in your data set specified in the argument data.It has a similar structure as for functions in R for standard statistical methods such as lm for linear models, but it is more complicated.Let y, g, s, cc, and c be respectively the response, group, strata, categorical covariable, and continuous covariable names in the data.frame,and "r" is the reference group to be provided, the formula can be specified as follows: These terms can be omitted except for the response and group if you do not need them.In this formula, the strata and group variables, and the continuous and categorical covariables are recognized by functions strt, grp, covar and catecovar, respectively.In other words, all terms on the right hand side should be specified for the role by using these functions.For example, if one used the input sanon(y ~grp(g) + cc) with cc having no role specified, then cc would be ignored in the analysis.Covariables are distinguished on the basis of continuous or categorical type.The categorical covariable is internally converted into a set of dummy variables.The reference level for the group variable and the categorical covariable can be specified in the argument ref within each of the functions grp and catecovar.More than two strata variables can be in the formula object, and they are internally transformed into a cross-classification.The matrix P in (10) or ( 11) for the weighted least squares estimation is specified in the argument P by the R matrix object name.
The response variable can contain missing values, which should be coded by NA, and it can be multivariate (repeatedly measured).The function sanon has five options for dealing with missing values for the response variable, with the default maintaining missing values and using the methods described in Section 2.1.It can be specified by no argument or res.na.action = "default", and it assumes the MCAR mechanism for missing data.The second option is specified in function sanon by res.na.action = "LOCF1"; and it uses the LOCF convention for the kernels of U-statistics in (2), as described in Section 2.3.The function sanon recognizes the order in the left side of ~as a measurement order, and an available baseline can be the first value in the order.The third option is specified in the function sanon by res.na.action = "LOCF2"; and it applies the LOCF convention for Y as described in Section 2.3.The fourth is specified in the function sanon by res.na.action = "replace".The fifth option for managing missing data is the complete cases analysis, in which patients with missing values are removed, and this is specified in function sanon by res.na.action = "remove".
Support functions for objects returned by sanon are also available.First of all there is a print method for the display for the output of the function sanon.

print(object)
The object is an object of class 'sanon', usually, a result of a call to sanon, i.e., it is obtained by object <-sanon(...).Explicitly calling the print method is equivalent to typing object in the R console.The output consists of the input, the sample size, the strata, the response levels, the design matrix, and the stratified adjusted Mann-Whitney estimates.
The reference group and adjusting factor names are also represented to help interpretation.
A coef method is also available.It extracts the resulting adjusted estimates b defined in (10) from the object which is the output of sanon.

coef(object)
This is also a part of the output of the print method, but one can obtain it as a vector object to be available for additional algebraic calculations in R.
The corresponding estimated covariance matrix V b can be extracted using Because V b is not represented as the output of the print method, the vcov method can be helpful to check it.This returns a matrix object which is then available for additional algebraic calculations in R.
The results of inference after weighted least squares estimation can be obtained through the summary method.

summary(object)
The output is a p × 4 matrix with columns for the estimated coefficient, its standard error, chi-squared statistic and corresponding (two-sided) p value.The estimate in this function is ξ − 0.5 and the corresponding p value is for the hypothesis H 0 : ξ = 0.5.
A confint method for computing the confidence interval for b is also available.
confint(object, parm = NULL, level = 0.95) A specification of which parameters are to be given confidence intervals, either a vector of numbers or a vector of names, is needed.If missing, all parameters are considered.The confidence level can be specified.
The linear statistics c b for assessments of homogeneity or average across the response variables and the corresponding confidence intervals for the average are computed by the function contrast where c is the contrast matrix and b is a weighted least squares estimator in (10).This is used for this purpose after subtraction of 0.5 from c b.
This function provides the inference based on contrasts after applying the function sanon.
The contrast matrix C should be defined by the user.If C has only one row, and the input is specified as confint = TRUE (default FALSE), the confidence interval for the estimator is produced.
More details for the usage of function sanon will be explained through examples in Section 5.An alternative use -instead of the formula interface -is to specify R objects for arguments outcome, group, and strata variables and covariables; see help("sanon") for details.

Illustrative data sets
There are five data sets in the sanon package.The first two example data sets were used in Kawaguchi et al. (2011), and the other three example data sets were used in other papers.
As mentioned in the previous section, the following code works after installing the package using install.packages("sanon").

Chronic pain data
The first example data set in the sanon package is cpain, which can be considered as follows.
R> library("sanon") R> data("cpain", package = "sanon") R> summary(cpain) These data are from a multi-center randomized clinical trial to compare test and control treatments for the management of chronic pain, and they have been previously considered in Stokes et al. (2012, Chapter 14).This clinical trial has 8 strata (in correspondence to 2 centers × 4 diagnoses) for which the range of sample sizes is 10 to 34 and a univariate ordinal response variable with 5 categories (as excellent, good, moderate, fair, poor) for pain status after treatment for 4 weeks.

Respiratory disorder data
The second example data set in the sanon package is resp, which can be considered as follows.
R>  (1990).This clinical trial has 111 patients from two centers, and it has four post-baseline visits with corresponding ordinal response variables for patient global ratings of symptom control according to 5 categories (as 4 = excellent, 3 = good, 2 = fair, 1 = poor, 0 = terrible).

Relief of heartburn data
The third example data set in the sanon package is heartburn, which can be considered as follows.The data are from a two period crossover design clinical trial for relief of heartburn, and listings of the data appear in Koch, Gitomer, Skalland, and Stokes (1983).This clinical trial has 30 patients at each of two centers, 15 randomly assigned to the A:P sequence group (active treatment for the first period and placebo for the second period) and 15 to the P:A sequence group (placebo for the first period and active for the second period).The ordinal response (MD1 and MD2) is a composite measure for time to relief.This composite measure has values between 1 and 20 if relief occurred within 20 minutes from first dose of treatment; between 21 and 40 if relief occurred within 20 minutes of the second dose (and there was no relief from the first dose within 15 minutes and use of the second dose shortly thereafter); and the value 60 if there was no relief from either the first or second dose within the respective 20 minute limits.Koch et al. (1983) used the binary response which is whether or not the patient experienced relief within 15 minutes of the first dose of treatment (res1 and res2).

Seborrheic dermatitis data
The fourth example data set in the sanon package is sebor, which can be considered as follows.
R> data("sebor", package = "sanon") R> summary(sebor) R> head(sebor) The data are from a randomized clinical trial to compare a test treatment to placebo for seborrheic dermatitis, and listings of the data appear in Ramaswamy, Koch, and Amara (1997).This clinical trial has 167 patients from eight centers, and it has three anatomical sites (face, scalp, and chest) corresponding to ordinal response variables for patient global scores according to 6 categories (as 0 = cleared, 1 = excellent improvement, 2 = moderate improvement, 3 = slight improvement, 4 = no change, 5 = exacerbation).

Skin conditions data
The fifth example data set in the sanon package is skin, which can be considered as follows.
A confint method for computing the confidence interval is also available.
For illustrative purposes, diagnosis is alternatively managed as a categorical covariable with center remaining as a factor for stratification.This structure is of exploratory interest since it has substantially larger sample sizes within its strata.

M-W Estimate and 95% Confidence Intervals (adjusted by diagnosis[C/D], diagnosis[B/D], diagnosis[A/D] ):
Estimate Lower Upper response 0.5729 0.4971 0.6488 The result is very similar to the previous analysis with the adjusted Mann-Whitney estimator being 0.5729 (95%CI: 0.4971-0.6488;p = 0.059).As noted in Section 2.2, the management of a baseline factor as strata or as a covariable should be pre-specified.

Respiratory disorder data
The second example is a randomized clinical trial with four strata, two covariables, and four ordinal response variables.This example has gender (female or male) as an additional factor for stratification, and so there are 4 strata for center × gender; and it has age and the baseline rating of symptom control (with the same ordinal categories as the response variables) as two covariables.The application of function sanon to this data set is as follows.
The resulting vector of stratification adjusted estimators for the comparisons between the test treatment and placebo is in (13) containing ξ0 , ξ1 , ξ2 , ξ3 , ξ4 for the Mann-Whitney estimators that correspond to the baseline visit and visits 1, 2, 3, 4 and g for the difference between mean ages..4799, 0.6005, 0.7139, 0.6535, 0.6155, 1.0501] . ( The corresponding estimated covariance matrix from the methods in Section 2 can be obtained as follows. Both ( ξ0 − 0.5) and g have expected values of null on the basis of randomization of patients to the two treatment groups.The function contrast is used for this purpose after subtraction of 0.5 from ξ0 , ξ1 , ξ2 , ξ3 and ξ4 .
The resulting adjusted estimators b from the methods in Section 2.2 can be computed as follows.

R> summary(out22)
Call: sanon.formula(formula = cbind(baseline, visit1, visit2, visit3, visit4)  The respective p values are 0.0152, < 0.0001, 0.0012, 0.0155 for visits 1, 2, 3, and 4. Homogeneity of the ξ k across the four visits can be assessed with Since Q C 1 = 8.93 has p = 0.03 with respect to its approximate chi-squared distribution with d.f.= 3, there is some suggestion of departures of the ξ k from homogeneity; and from the inspection of their estimates, it appears that the difference between test treatment and placebo tends to be larger at visit 2 and visit 3 than at visit 1 and visit 4.
A comparison between treatments for the average of the ξk across the 4 visits is possible with Thus, the result for testing the null hypothesis H 0 : C 2 ξ = 0.5 is Q C 2 = 16 for which a two-sided p value of p < 0.0001 with respect to the approximate chi-squared distribution with d.f.= 1 is obtained.Note that contrast computes the Mann-Whitney estimate and its confidence interval (CI) adjusted for strata and covariables with the specification confint = TRUE in the case where C has one row (i.e., d.f.= 1).The estimate averaged over visits was C 2 b + 0.5 = 0.6548 (95%CI: 0.5789-0.7306).

Relief of heartburn data
The third example is a crossover design randomized clinical trial with two strata, one continuous covariable, and two ordinal response variables.The two strata are for centers, and age is a continuous covariable.For illustrative convenience, a missing value in age is imputed with the mean for the corresponding stratum.The statistics {b k /s.e.(b k )} 2 , where the b k are the adjusted estimators for (ξ k − 0.5) have approximately chi-squared distributions with d.f.= 1 under the null hypotheses H 0k : ξ k = 0.5.Although both b 1 for the comparison between P and A for the first period and b 2 for the comparison between P and A for the second period indicate more favorable responses for active treatment A than for placebo P, only the second period has p < 0.05 for the comparison between A and P. Since it is recognized that the comparison between the AP and PA groups for period 1 pertains to A versus P whereas that for period 2 pertains to P versus A, treatment×period interaction is assessed for
Since the two periods of this crossover study were separated by a sufficiently long washout period, the treatment × period interaction corresponding to this departure from homogeneity is unlikely to be due to unequal pharmacological carryover effects of A and P during the first period.
The Mann-Whitney estimates and their confidence intervals adjusted for strata and covariables for each period are computed as follows.

Seborrheic dermatitis data
The fourth example is a randomized clinical trial with seven strata (after the pooling of centers 6 and 7 in the original data as in Ramaswamy et al. 1997), three ordinal covariables with missing values, and three ordinal response variables with missing values.For this example, the MCAR assumption is realistic because missing sites did not receive treatment.
R> sebor2 <-sebor R> sebor2$center <-ifelse(sebor2$center == 6, 7, sebor2$center) This example has 7 centers and disease severities recorded at the baseline measurement.In the previous three examples, centers were managed as strata to account for the stratified randomization within centers.However, for this example they are managed as a categorical covariable, mainly because the sample size of 18 patients for score 3 is insufficient to support stratified estimation (relative to the 167 patients in this clinical trial).Additionally, initial disease severities can be considered as covariables; and they are ordinal with some missing values.Thus, they are treated in the left side of ~in the function sanon as responses.
To manage initial disease severities as covariables, the model P = [I 3 , 0 3,9 ] is fitted to f (after subtraction of 0.5 from each of the ξk ) by weighted least squares; here 0 3,9 is a 3×9 matrix with all entries being 0.
A comparison between treatments for the average of the ξk across the 3 sites is possible with the model P = [1 1 1 0 0 0 0 0 0 0 0 0] , and it is fitted to f (after subtraction of 0.5 from each of the ξk ) by weighted least squares.Its result for the overall comparison between treatments has a two-sided p value of p = 0.086.The estimate averaged over visits was b + 0.5 = 0.4338 (95%CI: 0.3581-0.5095).This comparison averages sites with more weight for those with larger sample size (and thereby less variance).

Skin conditions data
The fifth example is a randomized clinical trial with six strata, a categorical covariable with three levels, and three ordinal response variables with missing values.First of all, the centers 3 and 4 in the original data set are pooled, because there are only four observations in center 4, and rankings are less informative when overall sample size (n h 1 + n h 2 ) = n h within strata are < 10; see Kawaguchi et al. (2011).
R> skin2 <-skin R> skin2$center <-ifelse(skin2$center == 4, 3, skin2$center) This example has 5 strata for center and it has initial severity of the skin conditions as a categorical covariable.At three follow-up visits, patients were evaluated according to a fivepoint ordinal response scale defining extent of improvement.There are some missing values in response variables.The application of function sanon to this data set is as follows.Note that the estimates are for (ξ k − 0.5).Homogeneity of the ξ k across the three visits can be assessed with Contrast Inference: Chisq df Pr(>Chisq) [1,] 3.58 2 0.17 Since Q C 1 = 3.58 has p = 0.17 with respect to its approximate chi-squared distribution with d.f.= 2, there is no suggestion of departures of the ξk from homogeneity.
A comparison between treatments for the average of the ξk across the 3 visits is possible with Its result is Q C 2 = 152 for which a two-sided p value of p < 0.0001 with respect to the approximate chi-squared distribution with d.f.= 1 is obtained.The estimate averaged over visits was 0.1609 (95%CI: 0.1069-0.2149).In the case of opposite order of response variables, it was 0.8391 (95%CI: 0.7851-0.8931).
The function sanon has four other options for dealing with missing values for the response variable, although the default maintains missing values and uses the method described in Section 2 which can be selected by not specifying this argument or res.na.action = "default".
The second is the LOCF method based on the kernels of U-statistics in (2), which is specified in the function sanon by res.na.action = "LOCF1".In this example, the p value from method "LOCF1" was similar to that from the "default" method (i.e., the method described in Section 2).

R>
The third option is the LOCF method based on the observed value of Y , which is specified in the function sanon by res.na.action = "LOCF2".In this example, the p value from method "LOCF2" was similar to those from both "LOCF1" and "default".
The fourth method manages missing values as tied with all other values in the same stratum, which is specified in function sanon by res.na.action = "replace".The result shows that the p value after managing missing values as tied with all other values in the same stratum was similar to those for "LOCF1", "LOCF2", and "default".
The fifth option for managing missing data is the complete cases analysis, in which patients with missing values are removed, and it is specified in the function sanon by res.na.action = "remove".A comparison between treatments for the average of the ξk across the 3 visits is obtained using the following code.The result shows that the p value after removing subjects with missing values was similar to that for the method described in Section 2.
Note that the complete case analysis makes the same MCAR assumption as the default method, but has the limitation of only using complete cases rather than all of the data with missing maintained as missing.

Summary
The R package sanon contains functions to implement the methods for stratified Mann-Whitney estimators.These methods address the same comparisons between two randomized groups for a strictly ordinal response variable as the van Elteren test statistic.Since the functions have similar structures to standard R functions, they should be easily accessible for R users.The role of variables in the analysis can be specified in the formula by using functions from package sanon.For example, the stratification variable is specified by the function strt.Among these, functions treat and catecovar have arguments to set the reference group.The function sanon can deal with missing values in five ways, according to the user's intent.Function sanon is also applicable to data without stratification.The output has an orderly nature for the interpretation of the results.Thus, package sanon would be helpful for the analysis of randomized clinical trials with ordinal response variables.The details of the functions can be seen in the help files as listed using the command help(package = "sanon") in R.
(formula = response ~grp(treat, ref = "placebo") + strt(center) + strt(diagnosis), data = cpain) Sample size: 193 Strata ( center*diagnosis ): I*A, I*B, I*C, I*D, II*A, II*B, II*C, II*D Response levels: [response; 5 levels] (lower) poor, fair, moderate, good, excellent (higher)The argument data = cpain is for specification of the data set.Each variable in the data set is given a role by the formula which has a similar nature as the standard R function such as lm for linear models.The outcome response is put in the left side of ~.In the right side, the group variable treat and the strata variables center and diagnosis are specified by using the functions grp and strt, respectively, which are connected with +.The reference group of treat can be specified in the function grp by ref = "placebo".Two strata variables are taken as the cross-classification of two centers (as I, II) and four diagnoses (as A, B, C, D) corresponding to 8 strata.The result is stored in the object out11 and shown through the print method for 'sanon' objects, which indicates the resulting stratified Mann-Whitney estimate of ξ = 0.5804.The summary method produces the inference of the stratified Mann-Whitney estimator based on the null hypothesis of H 0 : ξ = 0.5.R> summary(out11)Call: sanon.formula(formula= response ~grp(treat, ref = "placebo") + strt(center) + strt(diagnosis), data = cpain) estimates of responses are for the (MW estimate -0.5).
R> contrast(out55, C = matrix(rep(1, 3)/3, ncol = 3), confint = TRUE) This example is from a randomized clinical trial to compare a test treatment to placebo for a respiratory disorder, and listings of the data appear inStokes et al. (2012, Chapter 15,  pp.515-516)andKoch,Carr, Amara, Stokes, and Uryniak The application of function sanon to this data set is as follows.