# Sensitivity Analysis for Principal Components of M-Scores for Several Outcomes in an Observational Study. For k outcomes in a matched observational study, principal() computes the M-scores for the outcomes one at a time, computes the principal components of the M-scores, and uses some of the larger principal components as the outcomes in a sensitivity analysis. The user controls the number of components using w: (i) if is.null(w)==TRUE, then the first principal component is used in a one-sided test, (ii) if length(w)==1, then w=1 and w=-1 both use the first principal component, but direct attention to the upper or lower tails, respectively, (iii) if length(w)>1, then the first length(w) principal components are used with weights w; e.g., w=c(1,1) adds the first two principal components together. Setting Scheffe=TRUE with length(w)=2 permits the user to test every linear combination of the first two principal components -- that is, every with length(w)=2 -- while controlling the family-wise error rate. Every matched set contains one treated subject and one or more controls. ```r principal(y,z,mset,w=NULL,gamma=1,inner=0,trim=3,lambda=0.5, TonT=FALSE,apriori=FALSE,Scheffe=FALSE,detail=FALSE, cor=FALSE) ``` ## Arguments - `y`: A matrix of responses with no missing data. Different columns of y are different variables, and there are k=dim(y)[2] variables. If present, the column names of y are used to label output. - `z`: Treatment indicators, z=1 for treated, z=0 for control with length(z)==dim(y)[1]. - `mset`: Matched set indicators, 1, 2, ..., sum(z) with length(mset)==dim(y)[1]. The vector mset may contain integers or may be a factor. - `w`: A vector of weights to be applied to principal components. There are k=dim(y)[2] variables in y. (i) If is.null(w)=TRUE or if length(w)=1 with w!=0, then the test is applied to the first principal component of the nvars M-test scores, and no adjustment for multiple testing is needed. (ii) If length(w)>1, then w determines a comparison among the first length(w) principal components of the k M-scores. At least one weight must be nonzero. The meaning of the weights is affected by whether cor=FALSE (the default) or cor=TRUE. If Scheffe=TRUE, the dimensionality of the Scheffe correction is length(w), so w=c(1,1) is different from w=(1,1,0), because the latter implies the investigator may consider the third principal component in some comparison, even though w=(1,1,0) ignores it. - `gamma`: gamma is the sensitivity parameter $\Gamma$, where $\Gamma \ge 1$. Setting $\Gamma = 1$ is equivalent to assuming ignorable treatment assignment given the matched sets, and it performs a within-set randomization test. - `inner`: inner and trim together define the $\psi$-function for the M-statistic. The default values yield a version of Huber's $\psi$-function, while setting inner = 0 and trim = Inf uses the mean within each matched set. The $\psi$-function is an odd function, so $\psi(w) = -\psi(-w)$. For $w \ge 0$, the $\psi$-function is $\psi(w)=0$ for $0 \le w \le $ inner, is $\psi(w)= $ trim for $w \ge $ trim, and rises linearly from 0 to trim for inner < w < trim. An error will result unless $0 \le $ inner $ \le $ trim. - `trim`: inner and trim together define the $\psi$-function for the M-statistic. See inner. Unlike other functions in this package, principal() requires trim1. In either case, this is a Normal approximation based on the central limit theorem and equals 1-pnorm(deviate). - `Scheffe`: If Scheffe=TRUE, then the weights are assumed to have been chosen after looking at the data. In this case, the P-value or P-value bound is corrected using Scheffe projections. The approximate corrected P-value or P-value bound is 1-pchisq(max(0,deviate)^2,length(w)). If Scheffe=FALSE and apriori=FALSE, then the deviate is returned, but no P-value is given. See Rosenbaum (2016). Note carefully that length(w) determines the extent of the correction for multiplicity, so that, for example, w=c(1,0) focuses attention on the first principal component but allows for consideration of all comparisons using the first two principal components. See the discussion of w above and the examples below. A Scheffe correction entitles you to look in both tails, which you do by considering both w and -w. See the planScheffe() function for a combination of an apriori and Scheffe comparisons. - `detail`: If detail=TRUE, then some detail from the princomp() function in the stats package is returned. - `cor`: If cor=FALSE, the principal components of the M-scores are computed from the covariance matrix, but if cor=TRUE, then they are computed from the correlation matrix. Because the columns of y were scaled using lambda and scored by the same $\psi$-function, they are on a common scale, and it is reasonable to compute the principal components from the covariance matrix with cor=FALSE, the default. Setting cor=TRUE standardizes the columns of y twice, so perhaps it is pointless. In any event, because the weights w are applied to the principal components, and the latter are affected by cor, it follows that the meaning of w is affected by the value of cor. ## Details If y has k columns for k outcomes, then comparison computes k M-scores, one for each outcome, computes principal components from these scores, combines the scores into a single comparison using w, and computes a one-sided, upper-tailed deviate for a randomization test or a sensitivity analysis, as described in Rosenbaum (2007, 2016). Outcomes are scaled using by the $\lambda$ quantile of the absolute differences before applying the $\psi$-function. In this sense, when trim1. - **ScheffePval**: If Scheffe=TRUE, the deviate is compared to the upper tail of the chi-square distribution to produce either a P-value for gamma=1 or an upper bound on the P-value for gamma>1. - **weights**: The weights are returned. ## References Huber, P. (1981) Robust Statistics. New York: John Wiley. (M-estimates based on M-statistics.) Maritz, J. S. (1979). A note on exact robust confidence intervals for location. Biometrika 66 163--166. (Introduces exact permutation tests based on M-statistics by redefining the scaling parameter.) Rosenbaum, P. R. (2007). Sensitivity analysis for m-estimates, tests and confidence intervals in matched observational studies. Biometrics 63 456-64. (R package sensitivitymv) Rosenbaum, P. R. (2013). Impact of multiple matched controls on design sensitivity in observational studies. Biometrics 69 118-127. (Introduces inner trimming, inner>0.) Rosenbaum, P. R. (2015). Two R packages for sensitivity analysis in observational studies. Observational Studies, v. 1. (Free on-line.) Rosenbaum, P. R. (2016) Using Scheffe projections for multiple outcomes in an observational study of smoking and periondontal disease. Annals of Applied Statistics, 10, 1447-1471. Rosenbaum, P. R. (2017) Combining planned and discovered comparisons in observational studies. Manuscript. Rosenbaum, P. R., & Rubin, D. B. (1985). The bias due to incomplete matching. Biometrics, 41, 103-116. ## Author(s) Paul R. Rosenbaum. ## Note For confidence intervals for individual outcomes, use function senmCI(). Under Fisher's hypothesis of no treatment effect, the principal components of the outcomes are unaffected by the treatment, so they may be used in randomization tests of no effect. However, this logic does not permit confidence intervals for the magnitude of effect on a principal component. The principal() function computes principal components of M-scores, not of the outcomes themselves. This has various implications. The M-scores share a common, resistant scaling, so it is reasonable to consider principal components of the covariance matrix of the M-scores. In contrast, M-tests computed from principal components of outcomes are not resistant to outliers because the components themselves are not resistant to outliers. M-scores add to zero within each matched set; see the example for the mscorev() function. In this specific and limited sense, variation among M-scores reflects variation within matched sets rather than variation between matched sets. For example, if the matched sets had been exactly matched for age, then the M-scores would be uncorrelated with age. In contrast, principal components of outcomes reflect both variation within and variation between matched sets. For instance, principal components of outcomes might be correlated with age even if the sets had been matched exactly for age. When the matched set size is variable, the M-scores incorporate variable weights, and the principal components are affected by these weights. For this reason, principal components of M-scores are more interpretable when every matched set has the same size, say matched pairs or matching 1-to-2, and they may be difficult to interpret if the set sizes vary widely, say 1-1 mixed with 1-5. In thinking about the relationship between outcomes and their M-scores, it can be helpful to examine the small, univariate example for the mscorev() function. ## Examples ```r # Please READ the documentation for artcog, and in particular # the distinction between simulated and actual data. # The dontrun section refers to the acutal data and # reproduces results in Rosenbaum (2017). # The example immediately below uses the simulated data, # and is simply a numerical illustration. data(artcog) attach(artcog) # Randomization test using the first principal component of the simulated data. principal(cbind(words,wordsdelay,animals),arthritis,mset,w=1,apriori=TRUE,detail=TRUE) # Randomization test exploring a contrast of the first two principal components. principal(cbind(words,wordsdelay,animals),arthritis,mset,w=c(1,-.1),Scheffe=TRUE) # Sensitivity analysis using the first principal component of the simulated data. principal(cbind(words,wordsdelay,animals),arthritis,mset,w=1,gamma=1.2,apriori=TRUE) amplify(1.2,c(1.5,2)) ## Not run: # For this illustration, obtain the actual data, # as described in the documentation for artcog. # An illustration from Rosenbaum (2017) follows. data(artcog) attach(artcog) # A randomization test using the first principal component for the three memory scores. # The loadings show that the first component gives positive weight to each memory score. principal(cbind(words,wordsdelay,animals),arthritis,mset,w=1,apriori=TRUE,detail=TRUE) # # The comparison above is insensitive to a bias of gamma=1.45 principal(cbind(words,wordsdelay,animals),arthritis,mset,w=1,gamma=1.45,apriori=TRUE,detail=TRUE) # # gamma=1.45 is an unobserved covariate that more than triples the odd of a poor memory score # and more than doubles the odds of arthritis. amplify(1.45,c(2,3,4)) # # Although the first principal component is insensitive to a bias of gamma=1.45, each # of the three individual variables is sensitive to a bias of gamma=1.45 senm(words,arthritis,mset,gamma=1.45) senm(wordsdelay,arthritis,mset,gamma=1.45) senm(animals,arthritis,mset,gamma=1.45) # # Although not particularly useful or enlightening in this one example, we can # explore all weighted combinations of the first two principle components, # correcting for multiple testing using Scheffe projections for dimension 2. # This would be more interesting in an example with 50 outcomes, where we # might want to reduce the dimensionality to 2 or 3 from 50, rather than to 1. # We will do calculations for gamma=1.25. A gamma=1.25 is an unobserved # covariate that doubles the odds of arthritis and doubles # the odds of a worse memory score. amplify(1.25,2) # The deviate is the same but the corrected P-value is different if w=1 or w=c(1,0), # because the former is doing a single one-sided test, while the latter is anticipating # consideration of all possible combinations of the first two components. principal(cbind(words,wordsdelay,animals),arthritis,mset,w=1,gamma=1.25,apriori=TRUE,detail=TRUE) principal(cbind(words,wordsdelay,animals),arthritis,mset,w=c(1,0),gamma=1.25,Scheffe =TRUE) # A weighted combination of the first two principal components, w=c(1,-.1), is ever so # slightly less sensitive than using the first component alone. principal(cbind(words,wordsdelay,animals),arthritis,mset,w=c(1,-.1),gamma=1.25,Scheffe =TRUE) detach(artcog) ## End(Not run) ```

principal function