How to evaluate the heat of formation with computational chemistry?
Your approach is quite correct, but as Jan already pointed out, it is incomplete. What you calculated is the difference in the electronic energy of the reaction $$\ce{H2 (g) + 1/2 O2 (g,\,{}^1\Delta_{g}) -&gt; H2O (g)}.$$

Let's turn this into some kind of a tutorial and a little exercise, that you can try to reproduce at home. For a more detailed description on how it works in Gaussian, you can read the paper by Joseph W. Ochterski. I will use Gaussian09 rev. D for the calculations involved.

The enthalpy of formation of water from the standard state of the elements involved, is basically the reaction enthalpy of the combustion of hydrogen and oxygen. You are in luck, because most of the quantum chemical programs shall provide you with the necessary data for reactions straight away. What you do need to do is calculate the thermal corrections to the reaction. Essentially this is done by calculating the Hessian, i.e. a frequency calculation. The second derivatives only have a interpretable meaning if they are performed at a stationary point. Hence they must be performed with the same method, that you used to optimise the geometry.
Jan already pointed out, that HF/6-31G(d) is probably not a well enough fit for geometries and frequencies. I'd always consider some level of theory with at least some coverage of correlation energy. I often prefer BP86 or M06L and the def2-SVP basis set as a first anchor point. The smaller the system, the more accurate the methods you can use (or afford).
Jan also pointed out, that especially for oxygen, the spin state is important. (I did a quick check and am almost 100% certain you calculated it in the singlet state; see above equation.)

For this exercise I have matched your suggestion with using M06/def2-QVZPP as it is still feasible. Keep in mind, that most of the times density functional approximations have to be calibrated. It's all about what you can or cannot afford. You should perform calculations to the best of your abilities.

Okay, let's dive into some mathematics. First, we consider the reaction equation 
$$\ce{H2 (g) + 1/2 O2 (g,\,{}^3\Sigma_{g}^{-}) -&gt; H2O (g)}.$$

And we would like to know the enthalpy of formation for water at standard state (which we will use as $298.15~\mathrm{K}$ and $1~\mathrm{atm}$)
$$\begin{align}
\Delta_{f}H^\circ(\ce{H2O}) &amp;= H^\circ(\ce{H2O}) 
                               - (H^\circ(\ce{H2}) 
                               +\frac12 H^\circ(\ce{O2}))\\
\Delta_{f}H^\circ(\ce{H2O}) &amp;= 
   [E_\mathrm{el}(\ce{H2O}) + H_\mathrm{corr}(\ce{H2O})]\\
&amp;\qquad
 - \left([E_\mathrm{el}(\ce{H2}) + H_\mathrm{corr}(\ce{H2})]
 +\frac12 [E_\mathrm{el}(\ce{O2}) + H_\mathrm{corr}(\ce{O2})]\right)\\
\end{align}$$

The data in the expanded form we can directly take from the output of the frequency calculation:

Oxygen  mo6qzvpp.freq.log/.com
 SCF Done:  E(UM06) =  -150.322374247     A.U. after    1 cycles
 Zero-point correction=                           0.003897 (Hartree/Particle)
 Thermal correction to Enthalpy=                  0.007203
 Thermal correction to Gibbs Free Energy=        -0.016049

Hydrogen mo6qzvpp.freq.log/.com
 SCF Done:  E(RM06) =  -1.17221845130     A.U. after    1 cycles
 Zero-point correction=                           0.009902 (Hartree/Particle)
 Thermal correction to Enthalpy=                  0.013206
 Thermal correction to Gibbs Free Energy=        -0.001590

Water mo6qzvpp.freq.log/.com
 SCF Done:  E(RM06) =  -76.4321360023     A.U. after    1 cycles
 Zero-point correction=                           0.021635 (Hartree/Particle)
 Thermal correction to Enthalpy=                  0.025415
 Thermal correction to Gibbs Free Energy=         0.004010


Therefore we will obtain 
$$\Delta_{f}H^\circ(\ce{H2O}) = -236.6~\mathrm{kJ/mol}$$
which I think is reasonably close to the experimental value.



Appendix
The following input files for the computation of above values.

%chk=mo6qzvpp.chk
#p M06/def2QZVPP
opt
int(ultrafinegrid)

oxygen

0 3
O
O 1 1.2
 (blank)




%chk=mo6qzvpp.freq.chk
%oldchk=mo6qzvpp.chk
#p M06/def2QZVPP
freq
geom=check guess=read
int(ultrafinegrid)

oxygen freq

0 3
 (blank)




%chk=mo6qzvpp.chk
#p M06/def2QZVPP
opt
int(ultrafinegrid)

hydrogen

0 1
H
H 1 1.2
 (blank)




%chk=mo6qzvpp.freq.chk
%oldchk=mo6qzvpp.chk
#p M06/def2QZVPP
freq
geom=check guess=read
int(ultrafinegrid)

hydrogen freq

0 1
 (blank)




%chk=mo6qzvpp.chk
#p M06/def2QZVPP
opt
int(ultrafinegrid)

water

0 1
O
H 1 1.2
H 1 1.2 2 109.0
 (blank)




%chk=mo6qzvpp.freq.chk
%oldchk=mo6qzvpp.chk
#p M06/def2QZVPP
freq
geom=check guess=read
int(ultrafinegrid)

water freq

0 1
 (blank)




Have a lot of fun.