  The extraction of cosmological parameters from microwave background
observations relies on specific assumptions about the statistical properties of
the data, in particular that the p-point distributions of temperature
fluctuations are jointly-normal. Using a battery of statistical tests, we
assess the multivariate Gaussian nature of the Wilkinson Microwave Anisotropy
Probe (WMAP) 1st year data. The statistics we use fall into three classes which
test different aspects of joint-normality: the first set assess the normality
of marginal (one-point) distributions using familiar univariate methods; the
second involves statistics that directly assess joint-normality; and the third
explores the evidence of non-linearity in the relationship between variates. We
applied these tests to frequency maps, `foreground-cleaned' assembly maps and
all-sky CMB-only maps. The assembly maps are of particular interest as when
combined with the kp2 mask, we recreate the region used in the computation of
the angular power spectrum. Significant departures from normality were found in
all the maps. In particular, the kurtosis coefficient, D'Agostino's statistic
and bivariate kurtosis calculated from temperature pairs extracted from all the
assembly maps were found to be non-normal at 99% confidence level. We found
that the results were unaffected by the size of the Galactic cut and were
evident on either hemisphere of the CMB sky. The latter suggests that the
non-Gaussianity is not simply related to previous claims of north-south
asymmetry or localized abnormalities detected through wavelet techniques.
