In genetic association studies, joint modeling of related traits (phenotypes) can utilize the correlation amongst them and thereby provide more power and uncover additional information about genetic etiology. Also, detecting rare genetic variants are of current scientific interest as a key to missing heritability. Logistic Bayesian LASSO (LBL) has been proposed recently to detect rare haplotype variants using case-control data, i.e., a single binary phenotype. As there is currently no haplotype association method that can handle multiple binary phenotypes, we extend LBL to fill this gap. We carry out extensive simulations to investigate the proposed bivariate LBL and compare with the univariate LBL. We find that the bivariate LBL performs better or similar to the univariate LBL in most settings. The bivariate LBL has the highest gain in power when there is a high correlation between the two traits and a haplotype associated with both traits affects at least one trait in a direction opposite to the direction of the correlation. We analyze Genetic Analysis Workshop 19 sequence data with phenotypes systolic and diastolic blood pressures and detect several associated rare haplotypes.