  Effects of possible orbital order in magnetic properties of two-dimensional
spin gap system for CaV$_4$O$_9$ are investigated theoretically. After
analyzing experimental data, we show that single orbital models assumed in the
literature are insufficient to reproduce the data. To understand the origin of
the discrepancy, we assume that in $d^1$ state of V, $d_{xz}$ and $d_{yz}$
orbitals have substantial contributions in the lowest-energy atomic level which
leads to a double-degeneracy. We study possible configurations of the orbital
order. By exact diagonalization and perturbation expansion, we calculate the
susceptibility, wavenumber dependence of low-lying excitations and equal-time
spin-spin correlations which is related to integrated intensity of the neutron
inelastic scattering. These quantities sensitively depend on the configuration
of the orbital order. The calculated results for some configurations of the
orbital order reproduce many experimental results much better than the previous
single-orbital models. However some discrepancy still remains to completely
reproduce all of the reported experimental results. To understand the origin of
these discrepancies, we point out the possible importance of the partially
occupied $d_{xy}$ orbital in addition to orbital order of partially filled
$d_{xz}$ and $d_{yz}$ orbitals.
