Neuroimaging phenotypes are often collected in genome-wide association studies (GWASs) as secondary phenotypes for a disease outcome. Joint analysis of multivariate imaging phenotypes can incorporate neural activity from multiple brain regions, and boost statis- tical power in association analysis by taking advantage of similarity across phenotypes. Yet, most GWASs are based on case-control study designs, implying that regression approaches not adjusted for the sampling design may lead to biased estimates of asso- ciations for secondary phenotypes with inflated Type I error rates and reduced power. Despite this well-known result, unadjusted regression models are widely used in the imaging genetic literature. The aim of this thesis is twofold: 1) to identify the condi- tions when sampling bias occurs in association analysis of secondary phenotypes, and 2) to improve power for gene discovery, utilizing multiple imaging phenotypes. Potential bias introduced by the unadjusted regression model is demonstrated using the Alzheimer’s Disease Neuroimaging Initiative (ADNI) data. In simulation studies, we compare the performance of the naive approach with those of several existing meth- ods accounting for ascertainment bias to demonstrate potential issues in using standard analyses of secondary phenotypes. Finally we propose two novel statistical methods to identify genetic associations with multiple phenotypes to improve testing power. The first method is to detect single-SNP and multi-trait associations in a proportional odds model (POM). The second considers multi-SNP and multi-trait associations in the gen- eralized estimating equations (GEE) framework, applied to rare variants in sequencing data and pathway analysis. Both methods extend the recently proposed adaptive sum of powered score (aSPU) test, shown to maintain high power in a wide range of situations. New methods are demonstrated in real data analyses and simulation studies.