Package 'sensitivityfull'

Title: Sensitivity Analysis for Full Matching in Observational Studies
Description: Sensitivity to unmeasured biases in an observational study that is a full match. Function senfm() performs tests and function senfmCI() creates confidence intervals. The method uses Huber's M-statistics, including least squares, and is described in Rosenbaum (2007, Biometrics) <DOI:10.1111/j.1541-0420.2006.00717.x>.
Authors: Paul R. Rosenbaum
Maintainer: Paul R. Rosenbaum <[email protected]>
License: GPL-2
Version: 1.5.6
Built: 2024-12-15 07:33:35 UTC
Source: CRAN

Help Index


Computes M-scores for a Full Match

Description

Of limited interest to most users, mscoref() computes the scores that form the basis for the hypothesis test performed by senfm.

Usage

mscoref(ymat, treated1, inner = 0, trim = 3, qu = 0.5)

Arguments

ymat

If there are I matched sets and the largest matched set contains J individuals, then y is an I by J matrix with one row for each matched set. If matched set i contains one treated individual and k controls, where k is at least 1 and at most J-1, then y[i,1] is the treated individual's response, y[i,2],...,y[i,k+1] are the responses of the k controls, and y[i,k+2],...,y[i,J] are equal to NA. If matched set i contains one control and k>1 treated individuals, then y[i,1] is the control's response, y[i,2],...,y[i,k+1] are the responses of the k treated individuals, and y[i,k+2],...,y[i,J] are equal to NA.

treated1

The vector treated1 is a logical vector of length I, where treated1[i]=TRUE if there is one treated subject in matched set i and treated1[i]=FALSE if there is more than one treated subject in matched set i.

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 ψ(w)=ψ(w)\psi(w) = -\psi(-w). For w >= 0, the ψ\psi-function is ψ(w)=0\psi(w)=0 for 0 <= w <= inner, is ψ(w)\psi(w) = trim for w >= trim, and rises linearly from 0 to trim for inner < w < trim.

trim

inner and trim together define the ψ\psi-function for the M-statistic. See inner.

qu

Before applying the ψ\psi-function to treated-minus-control differences, the differences are scaled by dividing by the qu quantile of all within set absolute differences. Typically, qu = 1/2 for the median. The value of qu has no effect if trim=0 and inner=Inf. See Maritz (1979) for the paired case and Rosenbaum (2007) for matched sets.

Value

Returns a matrix with the same dimensions as ymat and the same pattern of NAs. The returned value in position [i,j] compares ymat[i,j] to the other observations in row i of ymat, scoring the differences using ψ\psi-function, totalling them, and applying a weight. Matched sets of unequal size are weighted using weights that would be efficient in a randomization test under a simple model with additive set and treatment effects and errors with constant variance. See Rosenbaum (2007, section 4.2).

When a matched set contains one control and several treated subjects, this is reflected in the returned scores by a sign reversal.

Author(s)

Paul R. Rosenbaum

References

Huber, P. (1981) Robust Statistics. New York: John Wiley.

Maritz, J. S. (1979). A note on exact robust confidence intervals for location. Biometrika 66 163–166.

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. (2010). Design of Observational Studies. New York: Springer. Table 2.12, page 60, illustrates the calculations for the simple case of matched pairs.

Rosenbaum, P. R. (2013). Impact of multiple matched controls on design sensitivity in observational studies. Biometrics 69 118-127. (Introduces inner trimming to increase design sensitivity.)

Examples

# The artificial example that follows has I=9
# matched sets.  The first 3 sets have one treated
# individual and two controls with treated subjects
# in column 1.  The next 3 sets are
# matched pairs, with treated subjects in column 1.
# The next 3 sets have one control and two treated
# subjects, with the control in column 1.  Simulated
# from a Normal distribution with an additive effect
# of tau=1.

y<-c(2.2, 1.4, 1.6, 2.4, 0.7, 1.3, 1.2, 0.6, 0.3,
0.5, -0.1, -1.3, -0.3, 0.1, 0.4, 3.0, 1.1, 1.4, -0.8,
0.1, 0.8, NA, NA, NA, 1.1, 0.5, 1.8)
y<-matrix(y,9,3)
treated1<-c(rep(TRUE,6),rep(FALSE,3))

mscoref(y,treated1) # Huber scores
mscoref(y,treated1,inner=0.5,trim=3) #inner trimmed scores
mscoref(y,treated1,qu=.9,trim=1) #trimming the outer 10 percent

# For an additional example, install and load package sensitivitymv
# The following example is a match with variable controls.
# Both mscorev() (in sensitivitymv) and mscoref() (in sensitivityfull)
# reproduce the example in Rosenbaum (2007, Table 3).
# data(tbmetaphase)
# mscorev(tbmetaphase,trim=1)
# mscoref(tbmetaphase,rep(TRUE,15),trim=1)

Sensitivity Analysis for a Full Match in an Observational Study.

Description

In a full match, each matched set contains either one treated individual and one or more controls or one control and one or more treated individuals. Uses Huber's M-statistic as the basis for the test, for instance, a mean. Performs either a randomization test or an analysis of sensitivity to departures from random assignment. For confidence intervals, use function senfmCI().

Usage

senfm(y, treated1, gamma = 1, inner = 0, trim = 3, lambda = 1/2,
     tau = 0, alternative="greater")

Arguments

y

If there are I matched sets and the largest matched set contains J individuals, then y is an I by J matrix with one row for each matched set. If matched set i contains one treated individual and k controls, where k is at least 1 and at most J-1, then y[i,1] is the treated individual's response, y[i,2],...,y[i,k+1] are the responses of the k controls, and y[i,k+2],...,y[i,J] are equal to NA. If matched set i contains one control and k>1 treated individuals, then y[i,1] is the control's response, y[i,2],...,y[i,k+1] are the responses of the k treated individuals, and y[i,k+2],...,y[i,J] are equal to NA.

Although y can contain NA's, y[i,1] and y[i,2] must not be NA for every i.

treated1

The vector treated1 is a logical vector of length I, where treated1[i]=TRUE if there is one treated subject in matched set i and treated1[i]=FALSE if there is more than one treated subject in matched set i.

gamma

gamma is the sensitivity parameter Γ\Gamma, where Γ1\Gamma \ge 1. Setting Γ=1\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 ψ(w)=ψ(w)\psi(w) = -\psi(-w). For w0w \ge 0, the ψ\psi-function is ψ(w)=0\psi(w)=0 for 0w0 \le w \le inner, is ψ(w)=\psi(w)= trim for ww \ge trim, and rises linearly from 0 to trim for inner < w < trim.

An error will result unless 00 \le inner \le trim.

Taking trim < Inf limits the influence of outliers; see Huber (1981). Taking inner > 0 often increases design sensitivity; see Rosenbaum (2013).

trim

inner and trim together define the ψ\psi-function for the M-statistic. See inner.

lambda

Before applying the ψ\psi-function to treated-minus-control differences, the differences are scaled by dividing by the lambda quantile of all within set absolute differences. Typically, lambda = 1/2 for the median. The value of lambda has no effect if trim=0 and inner=Inf. See Maritz (1979) for the paired case and Rosenbaum (2007) for matched sets.

An error will result unless 0 < lambda < 1.

tau

The null hypothesis asserts that the treatment has an additive effect, tau. By default, tau=0, so by default the null hypothesis is Fisher's sharp null hypothesis of no treatment effect.

alternative

If alternative="greater", the null hypothesis of a treatment effect of tau is tested against the alternative of a treatment effect larger than tau. If alternative="less", the null hypothesis of a treatment effect of tau is tested against the alternative of a treatment effect smaller than tau. In particular, alternative="less" is equivalent to: (i) alternative="greater", (ii) y replaced by -y, and (iii) tau replaced by -tau. See the note for discussion of two-sided sensitivity analyses.

Details

For the given Γ\Gamma, senfm() computes the upper bound on the 1-sided P-value testing the null hypothesis of an additive treatment effect tau against the alternative hypothesis of a treatment effect larger than tau. By default, senfm() tests the null hypothesis of no treatment effect against the alternative of a positive treatment effect. The P-value is an approximate P-value based on a Normal approximation to the null distribution; see Rosenbaum (2007).

Matched sets of unequal size are weighted using weights that would be efficient in a randomization test under a simple model with additive set and treatment effects and errors with constant variance; see Rosenbaum (2007).

The upper bound on the P-value is based on the separable approximation described in Gastwirth, Krieger and Rosenbaum (2000); see also Rosenbaum (2007).

Value

pval

The upper bound on the 1-sided P-value.

deviate

The deviate that was compared to the Normal distribution to produce pval.

statistic

The value of the M-statistic.

expectation

The maximum expectation of the M-statistic for the given Γ\Gamma.

variance

The maximum variance of the M-statistic among treatment assignments that achieve the maximum expectation. Part of the separable approximation.

Note

The function senfm() performs 1-sided tests. One approach to a 2-sided, α\alpha-level test does both 1-sided tests at level α/2\alpha/2, and rejects the null hypothesis if either 1-sided test rejects. Equivalently, a bound on the two sided P-value is the smaller of 1 and twice the smaller of the two 1-sided P-values. This approach views a 2-sided test as two 1-sided tests with a Bonferroni correction; see Cox (1977, Section 4.2). In all cases, this approach is valid large sample test: a true null hypothesis is falsely rejected with probability at most α\alpha if the bias in treatment assignment is at most Γ\Gamma; so, this procedure is entirely safe to use. For a randomization test, Γ=1\Gamma=1, this Bonferroni procedure is not typically conservative. For large Γ\Gamma, this Bonferroni procedure tends to be somewhat conservative.

Related packages are sensitivitymv, sensitivitymw and sensitivity2x2xk. sensitivitymv is for matched sets with one treated subject and a variable number of controls. sensitivitymw is for matched sets with one treated subject and a fixed number of controls, including matched pairs. For their special cases, sensitivitymv and sensitivitymw contain additional features not available in sensitivityfull. sensitivitymw is faster and computes confidence intervals and point estimates. sensitivitymw also implements methods from Rosenbaum (2014). sensitivity2x2xk is for 2x2xk contingency tables, treatment x outcome x covariates; see Rosenbaum and Small (2016).

Rosenbaum (2007) describes the method for matching with variable numbers of controls, but only very minor adjustments are required for full matching, and senfm() implements these adjustments.

Author(s)

Paul R. Rosenbaum.

References

Cox, D. R. (1977). The role of signficance tests (with Discussion). Scand. J. Statist. 4, 49-70.

Hansen, B. B. (2007). Optmatch. R News 7 18-24. (R package optmatch) (Optmatch can create an optimal full match.)

Hansen, B. B. and Klopfer, S. O. (2006). Optimal full matching and related designs via network flows. J. Comput. Graph. Statist. 15 609-627. (R package optmatch)

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. (1991). A characterization of optimal designs for observational studies. J. Roy. Statist. Soc. B 53 597-610. (Introduces full matching.)

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

Rosenbaum, P. R. (2014). Weighted M-statistics with superior design sensitivity in matched observational studies with multiple controls. J. Am. Statist. Assoc. 109 1145-1158. (R package sensitivitymw)

Rosenbaum, P. R. and Small, D. S. (2016). An adaptive Mantel-Haenszel test for sensitivity analysis in observational studies. Biometrics, to appear.

Examples

# The artificial example that follows has I=9
# matched sets.  The first 3 sets have one treated
# individual and two controls with treated subjects
# in column 1.  The next 3 sets are
# matched pairs, with treated subjects in column 1.
# The next 3 sets have one control and two treated
# subjects, with the control in column 1.  Simulated
# from a Normal distribution with an additive effect
# of tau=1.

y<-c(2.2, 1.4, 1.6, 2.4, 0.7, 1.3, 1.2, 0.6, 0.3,
0.5, -0.1, -1.3, -0.3, 0.1, 0.4, 3.0, 1.1, 1.4, -0.8,
0.1, 0.8, NA, NA, NA, 1.1, 0.5, 1.8)
y<-matrix(y,9,3)
treated1<-c(rep(TRUE,6),rep(FALSE,3))

#Randomization test of no effect, Huber scores:
senfm(y,treated1)

#Sensitivity analysis, Huber scores:
senfm(y,treated1,gamma=2)

#Randomization test of tau=1 vs tau>1
senfm(y,treated1,tau=1)

#Randomization test of tau=1 vs tau<1
senfm(y,treated1,tau=1,alternative="less")

#Same randomization test of tau=1 vs tau<1
senfm(-y,treated1,tau=-1)

#Sensitivity analysis testing tau=1 at gamma=2
senfm(y,treated1,tau=1,gamma=2,alternative="greater")
senfm(y,treated1,tau=1,gamma=2,alternative="less")

# For an additional example, install and load package sensitivitymv
# The following example is a match with variable controls.
# So this example has one treated subject per matched set.
# Both mscorev (in sensitivitymv) and mscoref (in sensitivityfull)
# reproduce parts of the example in Rosenbaum (2007, Section 4).
# data(tbmetaphase)
# senmv(tbmetaphase,gamma=2,trim=1)
# senfm(tbmetaphase,rep(TRUE,15),gamma=2,trim=1)
# senmv(tbmetaphase,gamma=2,trim=1,tau=0.94)
# senfm(tbmetaphase,rep(TRUE,15),gamma=2,trim=1,tau=.94)
# senmv(tbmetaphase,gamma=2,trim=1,tau=0.945)
# senfm(tbmetaphase,rep(TRUE,15),gamma=2,trim=1,tau=.945)
# mscoref(tbmetaphase,rep(TRUE,15),trim=1)

Sensitivity Analysis for a Confidence Interval in a Full Match.

Description

In a full match, each matched set contains either one treated individual and one or more controls or one control and one or more treated individuals. Uses Huber's M-statistic as the basis for a confidence interval for an additive constant treatment effect, τ\tau. Performs either a randomization inference or an analysis of sensitivity to departures from random assignment. The confidence interval inverts the test in the function senfm() in the sensitivityfull package.

Usage

senfmCI(y,treated1,gamma=1,inner=0,trim=3,lambda=1/2,
                  alpha=0.05,twosided=TRUE,upper=TRUE)

Arguments

y

If there are I matched sets and the largest matched set contains J individuals, then y is an I by J matrix with one row for each matched set. If matched set i contains one treated individual and k controls, where k is at least 1 and at most J-1, then y[i,1] is the treated individual's response, y[i,2],...,y[i,k+1] are the responses of the k controls, and y[i,k+2],...,y[i,J] are equal to NA. If matched set i contains one control and k>1 treated individuals, then y[i,1] is the control's response, y[i,2],...,y[i,k+1] are the responses of the k treated individuals, and y[i,k+2],...,y[i,J] are equal to NA. Although y may, and typically does, contain NA's, y[i,1] and y[i,2] must not be NA for all i.

If you have matched pairs, not matched sets, use the senmwCI function in the sensitivitymw package.

treated1

The vector treated1 is a logical vector of length I, where treated1[i]=TRUE if there is one treated subject in matched set i and treated1[i]=FALSE if there is more than one treated subject in matched set i.

gamma

gamma is the sensitivity parameter Γ\Gamma, where Γ1\Gamma \ge 1. Setting Γ=1\Gamma = 1 is equivalent to assuming ignorable treatment assignment given the matched sets, and it returns a randomization-based confidence interval.

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 ψ(w)=ψ(w)\psi(w) = -\psi(-w). For w0w \ge 0, the ψ\psi-function is ψ(w)=0\psi(w)=0 for 0w0 \le w \le inner, is ψ(w)=\psi(w)= trim for ww \ge trim, and rises linearly from 0 to trim for inner < w < trim.

An error will result unless 00 \le inner \le trim.

Taking trim < Inf limits the influence of outliers; see Huber (1981). Taking inner > 0 often increases design sensitivity; see Rosenbaum (2013).

trim

inner and trim together define the ψ\psi-function for the M-statistic. See inner.

lambda

Before applying the ψ\psi-function to treated-minus-control differences, the differences are scaled by dividing by the lambda quantile of all within set absolute differences. Typically, lambda = 1/2 for the median. The value of lambda has no effect if trim=0 and inner=Inf. See Maritz (1979) for the paired case and Rosenbaum (2007) for matched sets.

An error will result unless 0 < lambda < 1.

alpha

The coverage rate of the confidence interval is 1-alpha. If the bias in treatment assignment is at most Γ\Gamma, then the confidence interval will cover the true τ\tau with probability at least 1α1-\alpha.

twosided

If twosided==TRUE, then a two-sided 1α1-\alpha confidence interval is constructed. If twosided==FALSE, then a one-sided 1α1-\alpha confidence interval is constructed. The two sided interval is the intersection of two one-sided 1α/21-\alpha/2 intervals.

upper

If twosided==TRUE, then upper is ignored. If twosided==FALSE and upper=TRUE, then an upper 1α1-\alpha confidence interval is constructed. If twosided==FALSE and upper=FALSE, then a lower 1α1-\alpha confidence interval is constructed.

Details

For the given Γ\Gamma, senfmCI() inverts the test in the function senfm() to produce the confidence interval. That is, it tests every τ\tau and retains the values not rejected at level α\alpha.

The test is a large sample approximation based on a Normal approximation to the null distribution; see Rosenbaum (2007).

Matched sets of unequal size are weighted using weights that would be efficient in a randomization test under a simple model with additive set and treatment effects and errors with constant variance; see Rosenbaum (2007).

The upper bound on the P-value is based on the separable approximation described in Gastwirth, Krieger and Rosenbaum (2000); see also Rosenbaum (2007).

Value

PointEstimates

The interval of point estimates of τ\tau. If gamma=1, then the interval is a single point estimate.

ConfidenceInterval

The confidence interval for τ\tau.

description

Reminder of the coverage rate, gamma, and type of interval.

Note

In a sensitivity analysis, a one-sided confidence interval is not conservative; however, two-sided intervals formed as the intersection of two one-sided 1α/21-\alpha/2 intervals are somewhat conservative. See the discussion of two-sided tests in the documentation for senfm().

Author(s)

Paul R. Rosenbaum.

References

Hansen, B. B. (2007). Optmatch. R News 7 18-24. (R package optmatch) (Optmatch can create an optimal full match.)

Hansen, B. B. and Klopfer, S. O. (2006). Optimal full matching and related designs via network flows. J. Comput. Graph. Statist. 15 609-627. (R package optmatch)

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. (1991). A characterization of optimal designs for observational studies. J. Roy. Statist. Soc. B 53 597-610. (Introduces full matching.)

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

Rosenbaum, P. R. (2014). Weighted M-statistics with superior design sensitivity in matched observational studies with multiple controls. J. Am. Statist. Assoc. 109 1145-1158. (R package sensitivitymw)

Examples

# The artificial example that follows has I=9
# matched sets.  The first 3 sets have one treated
# individual and two controls with treated subjects
# in column 1.  The next 3 sets are
# matched pairs, with treated subjects in column 1.
# The next 3 sets have one control and two treated
# subjects, with the control in column 1.  Simulated
# from a Normal distribution with an additive effect
# of tau=1.

y<-c(2.2, 1.4, 1.6, 2.4, 0.7, 1.3, 1.2, 0.6, 0.3,
0.5, -0.1, -1.3, -0.3, 0.1, 0.4, 3.0, 1.1, 1.4, -0.8,
0.1, 0.8, NA, NA, NA, 1.1, 0.5, 1.8)
y<-matrix(y,9,3)
treated1<-c(rep(TRUE,6),rep(FALSE,3))

# Randomization interval and point estimate, Huber scores:
senfmCI(y,treated1)

# Uses senfm() to show how senfmCI() inverts the test.
senfm(y,treated1,tau=0.6172307) #P-value is 0.025
senfm(y,treated1,tau=2.0612746,alternative = "less") #P-value is 0.025
senfm(y,treated1,tau=1.345436) #Statistic is 0

senfmCI(y,treated1,gamma=1.5) #Sensitivity of two-sided CI
# The next two calculations relate one and two-sided intervals
senfmCI(y,treated1,gamma=1.5,twosided=FALSE,upper=TRUE,alpha=0.025)
senfmCI(y,treated1,gamma=1.5,twosided=FALSE,upper=FALSE,alpha=0.025)

# If an estimator is approximately Normal, then +/- a standard
# error is approximately a 2/3 confidence interval.  Going the
# other way, people sometimes suggest looking at a 2/3
# confidence interval as analogous to +/- a standard error.
senfmCI(y,treated1,gamma=1.5,alpha=1/3)

# For an additional example, install and load package sensitivitymw
# library(sensitivitymw)
# The mercury data is 397 triples, 1 treated, 2 controls.
# It is the example in Rosenbaum (2014).
# data(mercury)
# help(mercury)
# In this balanced design, senmwCI() and senfmCI() give the same CI.
# senmwCI(mercury,gamma=3)
# senfmCI(mercury,rep(TRUE,397),gamma=3,twosided=FALSE)

Computes the Separable Approximation.

Description

Of limited interest to most users, separable1f() is called by the main function, senfm(), in computing the sensitivity analysis. separable1f() is given scores produced by mscoref() and computes the separable approximation to the upper bound on the P-value.

Usage

separable1f(ymat, gamma = 1)

Arguments

ymat

A matrix of scores produced by mscoref.

gamma

The sensitivity parameter Γ1\Gamma \ge 1.

Details

See Gastwirth, Krieger and Rosenbaum (2000) and Rosenbaum (2007, section 4) for discussion of the separable approximation.

Value

pval

The upper bound on the 1-sided P-value.

deviate

The deviate that was compared to the Normal distribution to produce pval.

statistic

The value of the score-statistic.

expectation

The maximum expectation of the score-statistic for the given Γ\Gamma.

variance

The maximum variance of the score-statistic among treatment assignments that achieve the maximum expectation. Part of the separable approximation.

Author(s)

Paul R. Rosenbaum

References

Gastwirth, J. L., Krieger, A. M. and Rosenbaum, P. R. (2000). Asymptotic separability in sensitivity analysis. J. Roy. Statist. Soc. B. 62 545-555.

Rosenbaum, P. R. (2007). Sensitivity analysis for m-estimates, tests and confidence intervals in matched observational studies. Biometrics 63 456-64. (See section 4.) (R package sensitivitymv)

Examples

# The artificial example that follows has I=9
# matched sets.  The first 3 sets have one treated
# individual and two controls with treated subjects
# in column 1.  The next 3 sets are
# matched pairs, with treated subjects in column 1.
# The next 3 sets have one control and two treated
# subjects, with the control in column 1.  Simulated
# from a Normal distribution with an additive effect
# of tau=1.

y<-c(2.2, 1.4, 1.6, 2.4, 0.7, 1.3, 1.2, 0.6, 0.3,
0.5, -0.1, -1.3, -0.3, 0.1, 0.4, 3.0, 1.1, 1.4, -0.8,
0.1, 0.8, NA, NA, NA, 1.1, 0.5, 1.8)
y<-matrix(y,9,3)
treated1<-c(rep(TRUE,6),rep(FALSE,3))

# The same calculation done in two equivalent ways.
separable1f(mscoref(y,treated1),gamma=2)
senfm(y,treated1,gamma=2)