We combine a smoothing spline analysis of variance (SS-ANOVA) model and a log-linear model to build a partly flexible model for multivariate Bernoulli data. The joint distribution conditioning on the predictor variables is estimated. The log odds ratio is used to measure the association between outcome variables. A numerical scheme based on the block one-step successive over relaxation SOR–Newton-Ralphson algorithm is proposed to obtain an approximate solution for the variational problem. We extend the generalized approximate cross validation (GACV) and the randomized GACV for choosing smoothing parameters to the case of multivariate Bernoulli responses. The randomized version is fast and stable to compute and is used to adaptively select smoothing parameters in each block onestep SOR iteration. Approximate Bayesian confidence intervals are obtained for the flexible estimates of the conditional logit functions. Simulation studies are conducted to check the performance of the proposed method, using the comparative Kullback–Leibler distance as a yardstick. Finally, the model is applied to two-eye observational data from the Beaver Dam Eye Study, to examine the association of pigmentary abnormalities and various covariates.