2 Data analysis
We will show how to use the CovRegRF package on a generated data set. The data set consists of two multivariate data sets: \(\mathbf{X}_{n \times 3}\) and \(\mathbf{Y}_{n \times 3}\). The sample size (\(n\)) is 200. The covariance matrix of \(\mathbf{Y}\) depends on \(X_1\) and \(X_2\) (i.e. \(X_3\) is a noise variable). We load the data and split it into train and test sets:
library(CovRegRF)
data(data)
xvar.names <- colnames(data$X)
yvar.names <- colnames(data$Y)
data1 <- data.frame(data$X, data$Y)
set.seed(4567)
smp <- sample(1:nrow(data1), size = round(nrow(data1)*0.6), replace = FALSE)
traindata <- data1[smp,, drop=FALSE]
testdata <- data1[-smp, xvar.names, drop=FALSE]Firstly, we check the global effect of \(\mathbf{X}\) on the covariance matrix estimates by applying the significance test for the three covariates.
formula <- as.formula(paste(paste(yvar.names, collapse="+"), ".", sep=" ~ "))
globalsig.obj <- significance.test(formula, traindata, params.rfsrc = list(ntree = 200),
nperm = 10, test.vars = NULL)
globalsig.obj$pvalue
#> [1] 0Using 10 permutations, the estimated \(p\)-value is 0 which is smaller than the significance level (\(\alpha\)) of 0.05 and we reject the null hypothesis indicating the conditional covariance matrices significantly vary with the set of covariates. When performing a permutation test to estimate a \(p\)-value, we need more than 10 permutations. Using 500 permutations, the estimated \(p\)-value is 0.012. The computational time increases with the number of permutations.
Next, we apply the proposed method with covregrf() and
get the out-of-bag (OOB) covariance matrix estimates for the training
observations.
covregrf.obj <- covregrf(formula, traindata, params.rfsrc = list(ntree = 200))
pred.oob <- covregrf.obj$predicted.oob
head(pred.oob, 2)
#> [[1]]
#> y1 y2 y3
#> y1 1.1286879 0.8387699 1.101836
#> y2 0.8387699 1.9878143 1.507416
#> y3 1.1018365 1.5074164 3.591839
#>
#> [[2]]
#> y1 y2 y3
#> y1 1.353311 1.050629 1.761897
#> y2 1.050629 2.153400 2.142313
#> y3 1.761897 2.142313 4.453531Then, we get the variable importance (VIMP) measures for the covariates. VIMP measures reflect the predictive power of \(\mathbf{X}\) on the estimated covariance matrices. Also, we can plot the VIMP measures.
vimp.obj <- vimp(covregrf.obj)
vimp.obj$importance
#> x1 x2 x3
#> 1.38012921 0.51122499 0.07159317
plot.vimp(vimp.obj)
From the VIMP measures, we see that \(X_3\) has smaller importance than \(X_1\) and \(X_2\). We apply the significance test to evaluate the effect of \(X_3\) on the covariance matrices while controlling for \(X_1\) and \(X_2\).
partialsig.obj <- significance.test(formula, traindata, params.rfsrc = list(ntree = 200),
nperm = 10, test.vars = "x3")
partialsig.obj$pvalue
#> [1] 0.3Using 10 permutations, the estimated p-values is 0.3 and we fail to reject the null hypothesis, indicating that we do not have enough evidence to prove that \(X_3\) has an effect on the estimated covariance matrices while \(X_1\) and \(X_2\) are in the model. Using 500 permutations, the estimated \(p\)-value is 0.218.
Finally, we can get the covariance matrix predictions for the test observations.
pred.obj <- predict(covregrf.obj, testdata)
pred <- pred.obj$predicted
head(pred, 2)
#> [[1]]
#> y1 y2 y3
#> y1 1.0710647 0.5039413 0.6859008
#> y2 0.5039413 1.6077176 1.1050398
#> y3 0.6859008 1.1050398 2.7710218
#>
#> [[2]]
#> y1 y2 y3
#> y1 1.294158 1.306257 1.854186
#> y2 1.306257 2.183658 2.424231
#> y3 1.854186 2.424231 4.387248