mvctm (version 1.2)

mvctm: Multivariate Variance Components Tests for Multilevel Data

Description

This function performs a permutation test for a variance component for 2-level, 3-level or 4-level data. The response can be univariate or multivariate.

Usage

mvctm(fixed, cluster, data, leveltested, method = "ls", npermut = 1000, 
weight = "observation", affequiv = TRUE)

Arguments

fixed

An object of class ``formula" describing the fixed effects part of the model using the variables in the data frame data.

cluster

A vector giving the name of the variables in the data frame data to specify the clustering configuration. The order is important. For 2-level data it is a vector of dimension 1 specifying the level 1 cluster. For 3-level data, it is a vector of dimension 2. The first element specifies the level 1 (outer) cluster and the second one specifies the level 2 (inner) cluster. For 4-level data, it is a vector of dimension 3. The first element specifies the level 1 (outer) cluster, the second one specifies the level 2 (middle) cluster, and the last one specifies the level 3 (inner) cluster.

data

A data frame containing the data.

leveltested

An integer giving the level to be tested. It must be 1 for 2-level data, 1 or 2 for 3-level data, and 1, 2 or 3 for 4-level data. It corresponds to the element in cluster.

method

The scores to be used. The four choices "ls", "mixed", "rank" and "sign" are available. The default is "ls". The choice "mixed" is only available for a univariate response.

npermut

The number of random permutation used to perform the test. The default is 1000.

weight

The weight function to be used. The three choices "pair", "observation" and "cluster" are available. The default is "observation".

affequiv

Whether or not we want to use the affine-equivariant version of the tests. This is only relevant for a multivariate response and method="rank" or "sign". The default is TRUE.

Value

A list with the following two elements:

pvalue

The p-value of the test.

statistic

The value of the test statistic computed on the original data.

Details

With method="ls", the fixed effects are estimated by ordinary least-squares. Then the test is performed on the residuals from this fit. With method="mixed", the fixed effects are estimated with a linear mixed model. Then the test is performed on the marginal (population) residuals from this fit. With method="rank", a rank-based method is used to estimate the fixed effects. Then the test is performed on the ranks of the residuals from this fit. Finally, with method="sign", a sign-based method is used to estimate the fixed effects. Then the test is performed on the signs of the residuals from this fit. For multivariate data, spatial ranks and signs are used.

With a univariate response, method="sign" is not recommended because the test might be liberal.

With, weight="pair", observations in larger clusters at the level leveltested will have more weights. With, weight="cluster", the same weight is given to each cluster at the level leveltested. As a compromise between these two, the default weight="observation" gives an equal weight to each individual observation, with respect to the clusters at level leveltested.

References

Larocque, D., Nevalainen, J. and Oja, H. (2018). Multivariate Variance Components Tests for Multilevel Data. Les Cahiers du GERAD G-2018-58.

Examples

# NOT RUN {
data(toydata)

# Bivariate 2-level model.
# Classroom as the clusters. 
# Only an intercept is in the fixed part of the model.
# Test based on 200 permutations
mvctm(fixed=cbind(y1,y2)~1,cluster=c("classroom"),
data=toydata,leveltested=1,npermut=200)

# Same as above but The two covariates are in the fixed part of the model.
# Test based on 1000 permutations (default).
# }
# NOT RUN {
mvctm(fixed=cbind(y1,y2)~x1+x2,cluster=c("classroom"),
data=toydata,leveltested=1)
# }
# NOT RUN {
# Same as above but the rank scores are used.
# }
# NOT RUN {
mvctm(fixed=cbind(y1,y2)~x1+x2,cluster=c("classroom"),
data=toydata,leveltested=1, method="rank")
# }
# NOT RUN {
# Univariate 4-level model. 
# Classrooms, nested within schools, nested within regions.
# The variance component at the region level is tested. 
# The fixed effects are estimated with a linear mixed model.
# }
# NOT RUN {
mvctm(fixed=y1~x1+x2,cluster=c("region","school","classroom"),
data=toydata,leveltested=1,method="mixed")
# }
# NOT RUN {
# Same as above but the variance component at the school level is tested.
# }
# NOT RUN {
mvctm(fixed=y1~x1+x2,cluster=c("region","school","classroom"),
data=toydata,leveltested=2,method="mixed")
# }
# NOT RUN {
# Same as above but the variance component at the classroom level is tested.
# }
# NOT RUN {
mvctm(fixed=y1~x1+x2,cluster=c("region","school","classroom"),
data=toydata,leveltested=3,method="mixed")
# }
# NOT RUN {
# Univariate 3-level model.
# The variance component at the classroom level is tested.
# The fixed effects are removed with an M-estimator with the rlm function 
#      in the MASS package. 
# Then the residuals from this fit are used to perform the test. 
# The ~0 in the formula tells mvctm to use  mresid directly to perform 
#      the test without any centering or transformation.  
# }
# NOT RUN {
library("MASS")
toydata[,"mresid"]=rlm(y1~x1+x2,data=toydata)$residuals
mvctm(fixed=mresid~0,cluster=c("school","classroom"),
data=toydata,leveltested=2)
# }
# NOT RUN {
# }