For region-based sequencing data, power to detect genetic associations can be improved through analysis of multiple related phenotypes. With this motivation, in this paper, we propose a novel test to detect association simultaneously between a set of rare variants, such as those obtained by sequencing in a small genomic region, and multiple continuous phenotypes. We allow arbitrary correlations among the phenotypes, and build on a linear mixed model by assuming the effects of the variants following a multivariate normal distribution with mean zero and a specific covariance matrix. A data-adaptive variance component test based on score type of statistics is derived, in order to account for unknown correlation parameter in the covariance matrix of the variants effects. Since our approach calculates the p-value analytically, the proposed test procedure is computationally efficient. Extensive simulations and an application to the UK10K project show that our proposed multivariate test is generally more powerful than univariate tests, especially when there are pleiotropic effects or highly correlated phenotypes.